{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Day 4 Solutions" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from scipy.optimize import curve_fit" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "# working with the classes\n", "# NFW --> delta_sigma\n", "class constants:\n", " \"\"\"Useful constants\"\"\"\n", " G = 4.301e-9 #km^2 Mpc M_sun^-1 s^-2 gravitational constant\n", " H0 = 100. #h kms-1 Mpc-1 hubble constant at present\n", " omg_m = 0.315 #omega_matter\n", "\n", "class halo(constants):\n", " \"\"\"Useful functions for weak lensing signal modelling\"\"\"\n", " \n", " def __init__(self, log_M200m, con_par):\n", " self.m_tot = 10**log_M200m # h-1 Msun\n", " self.c = con_par # concentration parameter\n", " \n", " self.rho_crt = 3 * self.H0**2 /(8.0*np.pi*self.G) # rho critical\n", " \n", " self.r_200 = (3 * self.m_tot /(4*np.pi*200*self.rho_crt*self.omg_m))**(1./3.) # Mpc h-1\n", " self.rho_0 = self.c**3 * self.m_tot/(4 * np.pi * self.r_200**3 *(np.log(1 + self.c) - self.c/(1 + self.c)))\n", " #print(\"log_M200m = %s h-1 M_sun, c = %s \\n\"%(np.log10(self.m_tot), self.c))\n", " \n", " def nfw(self,r):\n", " \"\"\"given r, this gives nfw profile as per the instantiated parameters\"\"\" \n", " r_s = self.r_200/self.c\n", " value = self.rho_0/((r/r_s)*(1+r/r_s)**2)\n", " return value\n", " \n", " def sigma_nfw(self,r):\n", " \"\"\"analytical projection of NFW\"\"\"\n", " r_s = self.r_200/self.c\n", " k = 2*r_s*self.rho_0\n", " sig = 0.0*r\n", " c = 0\n", " for i in r:\n", " x = i/r_s\n", " if x < 1:\n", " value = (1 - np.arccosh(1/x)/np.sqrt(1-x**2))/(x**2-1)\n", " elif x > 1:\n", " value = (1 - np.arccos(1/x)/np.sqrt(x**2-1))/(x**2-1)\n", " else:\n", " value = 1./3.\n", " sig[c] = value*k\n", " c=c+1\n", " return sig\n", " \n", " def avg_sigma_nfw(self,r): \n", " \"\"\"analytical average projected of NFW\"\"\"\n", " r_s = self.r_200/self.c\n", " k = 2*r_s*self.rho_0\n", " sig = 0.0*r\n", " c=0\n", " for i in r:\n", " x = i/r_s\n", " if x < 1:\n", " value = np.arccosh(1/x)/np.sqrt(1-x**2) + np.log(x/2.0)\n", " value = value*2.0/x**2\n", " elif x > 1:\n", " value = np.arccos(1/x)/np.sqrt(x**2-1) + np.log(x/2.0)\n", " value = value*2.0/x**2\n", " else:\n", " value = 2*(1-np.log(2))\n", " sig[c] = value*k\n", " c=c+1\n", " return sig\n", " \n", " def esd(self,r): \n", " \"\"\"ESD profile from analytical predictions\"\"\"\n", " val = self.avg_sigma_nfw(r) - self.sigma_nfw(r)\n", " return val\n", "\n", " \n" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "hp = halo(14,5)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "14.0 5 11.443311952479494 1.1093952351546374 14.880882611979143\n" ] } ], "source": [ "print(np.log10(hp.m_tot), hp.c, np.log10(hp.rho_crt), hp.r_200, np.log10(hp.rho_0))" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, '$\\\\rho(r)$')" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAELCAYAAAAlTtoUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkpklEQVR4nO3dd3RV9ZrG8e+bRqihBjABAoRLEVEg0kWwgooVEFAUBREVR8dx7ujcmbnNW8YugtJFLCCgXhHrtdJCCSBFEAhFEnrvLclv/iDMzWQCJiTn7HP2eT5rsZZnn3P2flls8+RXtznnEBERKY4orwsQEZHwodAQEZFiU2iIiEixKTRERKTYFBoiIlJsCg0RESm2GK8LCKSaNWu6lJQUr8sQEQkrS5Ys2eOcq1XUe74OjZSUFDIyMrwuQ0QkrJjZz+d6T91TIiJSbAoNEREpNl+Ghpn1MrOxBw8e9LoUERFf8WVoOOc+ds4NTUhI8LoUERFf8WVoiIhIYCg0RESk2BQa5/DJiu3sPXLS6zJEREKKr9dpXKhdh04wfMpSos24pnlt+l6eTNcmtYiJVsaKSGRTaBQhsUo8Xz7elWkZWXywdCuf/7iDxMrluKNtMn3aJtOoViWvSxQR8YT5+cl9aWlprrQrwk/n5vHNT7uYnpHFt2t3k5vnaJdSnT5pydxwSV0qllPuioi/mNkS51xake8pNIpv16ETvL90K9Mzsti45ygV46K5qdVF9L08mTb1q2FmZXYtERGvKDTKmHOOJT/vZ1pGFrNWbOfYqVxSEyvR7/J63NY6iRqVypX5NUVEgkWhEUBHTubw6YrtTF28haVbDhAbbVzbojZ3Xl6fLqk1iY5S60NEwktYhoaZNQJ+AyQ453rnH4sC/ghUATKcc2+e7xzBCI2C1u08zHuLs/hgaTb7j50mqWp5erdNpk9aMsnVKgStDhGR0giZ0DCzicBNwC7nXMsCx3sArwDRwHjn3F8LvDejQGjcBtwC7AM+cc59fb7rBTs0zjqZk8tXq3cxdfEW5mbuAeCKJrW4M60e17aoTVyMpu6KSOg6X2gEe+rPJGAkMPnsATOLBkYB1wLZwGIzm+mcW13E95sC6c65MWY2AzhvaHilXEw0N7aqy42t6pK17xjTl2QzPSOLR95dSo2KcfRum0y/dvVpWLOi16WKiJRIUEPDOTfbzFIKHW4HZDrnNgKY2VTOtCaKCo1s4FT+f+cWdQ0zGwoMBahfv34ZVF069apX4Ilrf8VjVzdh9vrdTF20hfFzNzFm9kY6NqpBv3b16NGyDuVior0uVUTkF4XCIoMkIKvA62ygvZnVAP4EtDazp51zfwE+AF41syuA2UWdzDk3FhgLZ7qnAlp5CURHGd2bJtK9aSK7Dp1g+pJspi7ewmNTf6BahVjuaHOm9ZGaqIWDIhK6QiE0ippe5Jxze4FhhQ4eAwb/4gnNegG9UlNTy6bCMpZYJZ5Huqfy0JWNmbdhD1MXZTFp/mbGz91Eu5Tq9G9fj54t6xIfq9aHiISWoM+eyu+emnV2INzMOgK/c85dn//6aYD8lkWpeDUQfiF2Hz7J+0uzmbpoC5v3HiOh/JnWx4D2an2ISHCFzOyp/GJS+L+hEQOsA64GtgKLgQHOuR9Le61wCo2z8vIcCzbt5d2FW/jixx2cznV0aFSdu9o34PqL62jmlYgEXMiEhplNAboBNYGdwG+dcxPM7AbgZc5MuZ3onPtTKa9ztnvqgfXr15euaA/tPnyS6UuyeHfhFrL3H6dmpTj6pNWj/+X1qV9D6z5EJDBCJjSCLRxbGkXJy3PMydzDOwt+5qs1O3GcWfdxV/v6XN0sUVu2i0iZirjQ8EtLoyjbDx7nvcVZTF2UxY5DJ6hTJZ5+7eoxoF19EqvEe12eiPhAxIXGWX5paRQlJ3/L9rcXbmH2ut3ERBnXX1yHgR0b0L5hde24KyIXLJRWhEsZiYmO4rqL63DdxXXYvOco7yz8mWkZ2XyycjtNEisxsGMDbmudROX4WK9LFREf8WVLw8/dU+dz4nQuM5dv4630n1m59SAV46K5rU0SAzuk0LROZa/LE5Ewoe6pCLQ86wCT03/m4xXbOJWTR7uG1RnYoQE9WtYhVgPnInIeCo0Itv/oKaZlZPH2wp/J2nec2lXKcXf7BvRvX5+aeliUiBRBoSHk5Tm+W7eLN+ZtZs76PcRFR3HTpXUZ1CmFVslVvS5PREJIxIVGpI5pFFfmriNMTt/M+0uyOXoql9b1qzKoUwo9W9bVinMRibzQOEstjfM7dOI0MzKymZy+mc17j1Gr8tmuq3okVtaaD5FIpdCQ88rLc3y/bjeT5m/m+3W7iY02erW6iPu7NKRlUoLX5YlIkGmdhpxXVJTRvVki3ZslsmH3ESbP38z0Jdl8sGwr7RtWZ3CXhlzdvDbRUVowKBLp1NKQIh08fpr3Fm/hzfk/s/XAcRrUqMB9nVLok1aPiuX0u4aIn0Vc95QGwstOTm4en/+4gwlzN7FsywEqx8fQ7/J63NspheRq2mlXxI8iLjTOUkujbC3dsp+Jczfx2aodAPS4uA6Dr2hIm/rVPK5MRMqSxjSkTLSpX402A6qx9cBxJqdvZsrCLXyycjtpDaoxtGsjrmlemyiNe4j4mloacsGOnsxhWkYWE+ZuInv/cRrVrMiQKxpxe5skPd9cJIype0oCKic3j09X7WDs7A2s2nqImpXiuLdjCnd3aEC1inFelyciJRRxoaGBcG8450jfuJexszfy3drdlI+Npm9aMkOuaES96ho0FwkXERcaZ6ml4Z21Ow4zdvZGZi7fSm6eo+cldRnWtTGXJGuxoEioU2iIZ3YcPMEb8zfx7oItHD6ZwxVNavJwt1Q6NNLTBUVClUJDPHfoxGneXvAzE+duYs+RU7SuX5VHuqVyVbNEzbgSCTEKDQkZJ07nMj0ji9Hfb2TrgeM0rV2Zh7o15qZWdYnRw6FEQoJCQ0LO6dw8Pl6+jde/28D6XUeoV708D3ZtTO+2yZquK+IxhYaErLw8x1drdjLquw0szzpArcrlGNKlIXd3aKA9rkQ8otCQkOecI33DXl77bgNzM/dQrUIsQ65oxD0dG1A5Ptbr8kQiikJDwsrSLfsZ8fV6vlu7m4TysQzp0pB7O6dQReEhEhQRFxpa3OcPy7MO8Oo36/lqzS4qx8dwf+eG3N+5IQkVFB4igRRxoXGWWhr+sGrrQUZ8vZ4vV++kcrkYBnVO4f7ODbVFiUiAKDTEF1ZvO8TIb9fz6codVIyL5t5OKQy5ohHVFR4iZUqhIb6ydsdhXv1mPZ+s3E7FuBju79KQwV0aklBe3VYiZUGhIb60fudhXvpqHZ+u3EFC+ViGdm3EoE4pmqorUkoKDfG1VVsP8tLf1/H1T7uoUTGOh7o15u4ODbRIUOQCKTQkIizdsp8Xv1zH3Mw9JFYux6NXpdL38nqUi1F4iJSEQkMiyoKNe3nhy7Us3ryfpKrleezqJtzeJkl7W4kU0/lCQ/8Xie90aFSDaQ92ZPL97ahZKY5fv7+Ca1+azacrt+PnX5JEgkGhIb5kZnT9VS3+9khnxg5sS2y08fA7S7n1tfmkb9jrdXkiYStkQ8PMGpnZBDObUeBYNzObY2ajzaybd9VJuDAzrru4Dp891pXnerdi16ET9B+3gEFvLGL1tkNelycSdoIaGmY20cx2mdmqQsd7mNlaM8s0s6cAnHMbnXODC53CAUeAeCA7OFWLH0RHGX3S6vHtk9349xuasWzLAW58dQ5PvPcDWfuOeV2eSNgI6kC4mXXlzA/9yc65lvnHooF1wLWcCYLFQH/n3Or892c453rn/3eUcy7PzGoDLzrn7jrf9TQQLudy8NhpXv9+A2/M24RzMLBjAx7pnqrV5SKE0EC4c242sK/Q4XZAZn7L4hQwFbjlHN/Py//P/UC5gBUqvpdQIZanejbj2ye7cWvri3hj3iaufPZbRn2byfFTuV6XJxKyQmFMIwnIKvA6G0gysxpmNhpobWZPA5jZ7WY2BngLGFnUycxsqJllmFnG7t27A127hLmLqpbn2d6X8vnjXWnfqAbPfbGWq174jg+XZZOXp5lWIoUFfZ2GmaUAswp0T/UBrnfODcl/PRBo55x7tLTXUveUlNSiTfv446zVrNx6kEuTE/jPm1qQllLd67JEgipkuqfOIRuoV+B1MrCtNCc0s15mNvbgwYOlKkwiT7uG1fnokc680OdSdhw6Qe/R6TzyzlINlovkC4WWRgxnBsKvBrZyZiB8gHPux9JeSy0NKY1jp3IYO3sjY77fSK5z3N+5IY90b6zHz4rvhUxLw8ymAOlAUzPLNrPBzrkcYDjwBbAGmFbawFBLQ8pChbgYHr/mV3z7ZDdualWX0d9voNtz3/HOwp/Jyc375ROI+JD2nhIpphXZB3hm1hoWbd5H09qV+c+bWtClSU2vyxIpcyHT0hAJZ62Sq/Legx14/a42HD+dy90TFjLsrSVk79d4h0QOX7Y0zKwX0Cs1NfWB9evXe12O+NDJnFzGz9nEyG8yyXOOR7qnMrRrIz3DQ3xBW6OLBMi2A8f506dr+GTFdupXr8B/3dSCq5snYmZelyZywdQ9JRIgF1Utz6gBbXh3SHvKxUQxZHIG909azKY9R70uTSQgfBkamj0lwdYptSafPnYF/3FjcxZv3s/1L83m2c9/4tipHK9LEylT6p4SKWO7Dp/gr5/9xAdLt1I3IZ7f3NicGy+pqy4rCRvqnhIJosTK8bzY9zJmDOtItQpxDH93GfdMXMRmdVmJDyg0RAIkLaU6Hz/ahd/ffDHLthzgupdnM/Kb9ZzK0cJACV++DA2NaUioiI4y7u2Uwtf/ciXXNq/N81+u44YRc1i0qfATAkTCg8Y0RILom5928p9/+5GtB45zZ1o9nurZjGp68JOEGI1piISIq5rV5u9PdOXBKxsxY2k2V7/4Pe8vycbPv7yJvyg0RIKsQlwMT/dszqxHu9CgRgX+Zfpy7hq/kI27j3hdmsgv8mVoaExDwkHzulV4f1gnnrm1JSu3HqTHy3N45SsNlEto05iGSAjYdfgEf/h4NbNWbKdZnco827sVrZKrel2WRCiNaYiEuMTK8Ywc0IZx96Sx7+gpbh01j79+9hMnTud6XZrI/6HQEAkh17aozd+fuJI+besx+vsN3PDKHDI2a3quhA6FhkiISSgfy3/3bsVbg9txMiePPmPS+d3MHzl6UvtYifcUGiIh6oomtfjyn7tyT4cGTJq/metfns28zD1elyURzpehodlT4hcVy8Xw+1taMu3BjsRGR3HX+IU89f4KDp047XVpEqE0e0okTJw4nctLX61j3OyNJFaO5y93XEL3polelyU+pNlTIj4QHxvN0z2b8+HDnalSPob73ljMv3+4UmMdElQKDZEwc2m9qswc3oUHuzZiyqIt9HxlDos1w0qCRKEhEobiY6N5+obmvDe0Iw5H3zHp/OWzNZzM0boOCSyFhkgYa9ewOp891pV+l9dnzPcbufnVefy4TRNAJHAUGiJhrlK5GP5y+yW8Mehy9h87xS0j5zHym/Xk5GoPKyl7Cg0Rn+jeLJEvHu9Kj5Z1eP7LdfQena6dc6XM+TI0tE5DIlW1inGMHNCGEf1bs2nPUW4YMYfJ6Zv1vA4pM1qnIeJTOw+d4NczVvD9ut1c1SyRZ3u3omalcl6XJWFA6zREIlDtKvFMuu9yfterBXMz99Dj5Tl8v26312VJmFNoiPiYmTGoc0M+eqQz1SrEcu/ERTwza7Wm5soFK3FomFlFM4sORDEiEhjN61bh40e7cE/HBoyfu4nbRs0nc5cGyaXkfjE0zCzKzAaY2Sdmtgv4CdhuZj+a2XNm1iTwZYpIacXHRvOHW1oy7p40th88zk2vzmHKoi0aJJcSKU5L41ugMfA0UMc5V885lwhcASwA/mpmdwewRhEpQ9e2qM3nj3clrUF1nv5gJcPeXsL+o6e8LkvCxC/OnjKzWCAXeMo59+dzfcY5F3J7NWv2lMi55eU5xs/dyHNfrKVGxXK8eOeldGpc0+uyJASUavaUc+60cy4PuOZ8nylFfSLigagoY2jXxnz4cGcqxEVz1/iFvPDlWnLz1F0l51aSgfBlZvZbM9OMKxEfaZmUwKx/6sIdbZJ59ZtMBoxbwM5DJ7wuS0JUSQKgHtAP2GZmH5nZH82sT4DqwswamdkEM5tR6HhFM1tiZjcF6toikaZCXAzP97mU5/tcyorsg9zwyhzmrNeaDvn/ih0azrm+zrnmQAPg90Am0L4kFzOziWa2y8xWFTrew8zWmlmmmT2Vf72NzrnBRZzm34BpJbmuiBRP77bJzBzemeoV47hn4iJeVHeVFFKcKbdW8LVz7qRzbqlz7k3n3JNFfeY8JgE9Cp0/GhgF9ARaAP3NrMU5arkGWA3sLOb1RKSEmtSuzEfDO9O7TTIj1F0lhRRryq2ZPWpm9QseNLM4M7vKzN4E7i3OxZxzs4HCjxhrB2TmtyxOAVOBW85xiu5AB2AA8IDGV0QCo0JcDM/1uZQX1F0lhRTnh24Pzky5nWJm281stZltAtYD/YGXnHOTSlFDEpBV4HU2kGRmNcxsNNDazJ4GcM79xjn3OPAuMC5/Vtf/YWZDzSzDzDJ279ZNLlIad+R3V9WodKa76oUv1+o5HRGuRLvc5q/ZqAMcds4duKALmqUAs5xzLfNf9wGud84NyX89EGjnnHv0Qs5fkNZpiJSN46dy+e3MVUzLyKZ9w+qM6N+a2lXivS5LAqRMdrk1s0eArZxZBf69mRU1SH0hsjkzM+usZGBbaU6o52mIlK3ycdE82/sf3VU3jphD+oa9XpclHijJmMCTQCvnXBJnuqy6mNnvyqCGxUATM2toZnGcmdY7szQndM597JwbmpCQUAblichZZ7urEsrHcveEhYydvUF7V0WYkoTGEWAXgHNuOzAYuL0kFzOzKUA60NTMss1ssHMuBxgOfAGsAaY5534syXmLuI5aGiIBcmZ2VReuv7g2f/70Jx5+ZylHTuZ4XZYESbHHNMzsYeBq4N+cc5n5YxNTnXMdAlhfqWhMQyRwnHOMn7OJv37+Eyk1KjBmYFtSEyt7XZaUgTIZ03DOvQa8A4w3s32cWdy31sz6aHt0kchjZjzQtRFvD27PweOnuXnkPGatKNVwpISBEq1zcM594JzrBiQCbYBvgE7AmLIv7cKpe0okeDo2rsGsR6+gWZ3KDH93Gc/MWs1pTcv1rRJNuQ036p4SCZ5TOXn8+dM1TJq/mXYNqzNyQGsSK2tabjgqk+4pEZHziYuJ4nc3X8zLd17GiuwD3DRiLhmbC28AIeHOl6Gh7ikR79zaOul/n9HRb+wCJqdv1rRcH/FlaGidhoi3mtetwkfDu3Dlr2rxXx/9yFPvr+RkTq7XZUkZ8GVoiIj3EsrHMu6eNIZ3T+W9jCz6j13ALu2WG/YUGiISMFFRxpPXN2XUgDas2X6YXiPn8kPWAa/LklLwZWhoTEMktNzYqi7vP9SJ2Ogo+o5J5/0l2V6XJBfIl6GhMQ2R0NPioirMHN6FtvWr8S/Tl/PHWau1zXoY8mVoiEhoql4xjsmD2zGoUwoT5m5i0BuLOXDslNdlSQkoNEQkqGKjz6zneLZ3KxZt2sfNI+exdsdhr8uSYvJlaGhMQyT09U2rx9QHO3DidC63vTaPz1ft8LokKQZfhobGNETCQ5v61fj40S78qnZlhr29hFHfZmohYIjzZWiISPioXSWeqUM7cOtlF/HcF2t5YtpyTpzWQsBQFeN1ASIi8bHRvHTnZaQmVuL5L9fx896jjBmYRq3K5bwuTQpRS0NEQoKZMfyqJrx2VxtWbz/EraPm8dOOQ16XJYUoNEQkpNxwSV2mP9iJnLw87nhtPl+v2el1SVKAL0NDs6dEwtslyQl89EgXGtWqxJDJGYybvVED5CHCl6Gh2VMi4a9OQjzTHuxIz5Z1+NOna3jq/ZWcytEKcq/5MjRExB/Kx0Uzsn8b/umqMzvlDpywkP1HtYLcSwoNEQlpUVHGE9c15ZV+l7Es6wC3vjaPzF1HvC4rYik0RCQs3HJZElOHduDoyRxuf20eCzbu9bqkiKTQEJGw0aZ+NT58uDOJVeIZOGEhHy7TFuvBptAQkbBSr3oF3h/WibQG1fnn95bzylfrNbMqiBQaIhJ2EirE8ub97bi9TRIvfbWOJ6ev0MyqIPHlNiJm1gvolZqa6nUpIhIgcTFRvNDnUhpUr8hLX61j24HjjB7YloTysV6X5mu+bGlonYZIZDAzHrumCS/2vZSMn/dxx+vzydp3zOuyfM2XoSEikeX2NslMvr89uw6d4LbX5vFD1gGvS/IthYaI+ELHxjX44OFOlI+Lpt/YdD3UKUAUGiLiG6mJlfnw4c40q1OFh95Zwvg5G70uyXcUGiLiKzUrlWPKAx24vkUdnvlkDX+ctZq8PE3JLSsKDRHxnfJx0Yy6qw2DOqUwYe4mHp26jJM5ehpgWfDllFsRkego47e9WnBR1Xj+/OlP7Dl8krH3pGlKbimppSEivmVmDO3amFf6XcbSLfvpM3o+2w4c97qssKbQEBHfu+WyJN68rx3bD5zg9tfm6zGypaDQEJGI0Cm1JtOGdcTh6PN6OvM37PG6pLAUsqFhZo3MbIKZzShwrLmZjTazGWb2kJf1iUj4aV63Ch883Jk6CfEMmriYmcu3eV1S2AlqaJjZRDPbZWarCh3vYWZrzSzTzJ4CcM5tdM4NLvg559wa59wwoC+QFrzKRcQvkqqWZ8awTlxWryr/NGWZ1nKUULBbGpOAHgUPmFk0MAroCbQA+ptZi3OdwMxuBuYCXweuTBHxs4QKsUwe3I4bLjmzluMPH2stR3EFNTScc7OBfYUOtwMy81sWp4CpwC3nOcdM51wn4K6i3jezoWaWYWYZu3fvLqvSRcRn4mPPPH/8vs4pTJy3icfe+0HbqxdDKIxpJAFZBV5nA0lmVsPMRgOtzexpADPrZmYjzGwM8GlRJ3POjXXOpTnn0mrVqhXw4kUkfEVFGf91Uwue6tmMj5dv4/5JizlyMsfrskJaKCzusyKOOefcXmBYoYPfAd/94gn1PA0RKSYzY9iVjalRMY6nPljJgHELeGPQ5dSoVM7r0kJSKLQ0soF6BV4nA6Wa0qDnaYhISfVJq8fYgW1Zt/MwvUen67kc5xAKobEYaGJmDc0sDugHzPS4JhGJQFc3r807Q9qz98hJ7nh9Pmu2axFgYcGecjsFSAeamlm2mQ12zuUAw4EvgDXANOfcj6W8Ti8zG3vw4MHSFy0iEaVtg+pMH9YJM+g7Jp1FmwrP3Yls5px/p5mlpaW5jIwMr8sQkTCUvf8Y90xcxNb9xxk5oA3XtqjtdUlBY2ZLnHNFroULhe6pMqeWhoiUVnK1CswY1olmdavw4FsZvLd4i9clhQRfhoYGwkWkLFSvGMe7Q9rTpUkt/u39lYz6NhM/984Uhy9DQ0SkrFQsF8P4e9K4+dKLeO6Ltfwhwp8EGArrNEREQlpcTBQv33kZNSrF8ca8zRw8fppn72hFTHTk/d7ty9DQ4j4RKWtnV49XrxDHC39fx6HjOYwc0Jr42GivSwsqX8akxjREJBDMjEevbsLvb76Yr9bs5L43FnP4xGmvywoqX4aGiEgg3dsphZfvvIxFm/dx1/iF7Dt6yuuSgsaXoaEptyISaLe2TmLswLas3XE4op497svQUPeUiATD1c1rM/n+duw6dJI+o9PZuPuI1yUFnC9DQ0QkWNo3qsGUoR04cTqXPqPTWbXV3z0cCg0RkVJqmZTAtGEdKRcTRf+xC3y9X5VCQ0SkDDSuVYkZD3WiVpVyDJywkG9+2ul1SQHhy9DQQLiIeOGiquWZ/mBHmtSuxNDJS5i5vFSPBgpJvgwNDYSLiFdqVCrHlAc60KZBNR6buoypi/y10aEvQ0NExEuV42N58752dG1Si6c+WMn4ORu9LqnMKDRERAKgfFw04+5Jo2fLOjzzyRpe+Wq9L3bIVWiIiARIXEwUr/ZvzR1tknnpq3X8+dM1YR8cvtywUEQkVMRER/Fc71ZULBfNuDmbOHIyl2dubUl0lHld2gXxZWhol1sRCSVRUcbvb76YSuVieO27DRw7lcPzfS4lNgy3Vg+/iotBs6dEJNSYGb/u0Yx/vb4pH/2wjYffWcqJ07lel1VivgwNEZFQ9Uj3VH5/88X8ffVOhryZwbFTOV6XVCIKDRGRILu3UwrP97mU+Rv2MHDCIg4eD59ncig0REQ80LttMiMHtGFF9gEGjFsQNs/kUGiIiHjkhkvqMvaeNDJ3HaHf2HR2HT7hdUm/SKEhIuKh7k0TeWPQ5WTtO86dYxaE/MOcFBoiIh7rlFqTtwa3Y8/hk/Qdk07WvmNel3ROvgwN7XIrIuEmLaU67zzQnsMnckL6KYC+DA2t0xCRcNQquSpTh3bgdG4efccsYO2Ow16X9P/4MjRERMJV87pVeO/BjkRHQb+xoff4WIWGiEiISU2sxLQHO1IhLob+4xawdMt+r0v6XwoNEZEQ1KBGRaYN60j1inEMHL+QBRv3el0SoNAQEQlZSVXLM+3BjtStWp5Bbyxi9rrdXpek0BARCWW1q8QzdWgHGtasxJA3M/hq9U5P61FoiIiEuJqVyjHlgfY0r1uZYW8v4fNV2z2rRaEhIhIGqlaI460h7WmVnMAj7y7j4+XbPKlDoSEiEiaqxMcyeXB72tavxmNTl/Hhsuyg1xDSoWFmjcxsgpnNKHDsVjMbZ2Yfmdl1XtYnIhJslcrFMOn+y+nQqAZPTFvOtMVZQb1+0EPDzCaa2S4zW1XoeA8zW2tmmWb2FIBzbqNzbnDBzznn/uacewAYBNwZtMJFREJEhbgYJg66nC6pNfn1+yt4Z+HPQbu2Fy2NSUCPggfMLBoYBfQEWgD9zazFL5znP/K/IyISceJjoxl3TxpXNUvkNx+uYtK8TUG5btBDwzk3G9hX6HA7IDO/ZXEKmArcUtT37Yz/Bj5zzi0NbLUiIqErPjaa0Xe35fqLa/O7j1czbvbGgF8zVMY0koCCHXPZQJKZ1TCz0UBrM3s6/71HgWuA3mY2rPCJzGyomWWYWcbu3d4vhBERCaS4mChGDmjDja3q8qdP1zDq28yAXi8moGcvPivimHPO7QWGFTo4AhhxrhM558YCYwHS0tJcWRYpIhKKYqOjeOXOy4iNMp77Yi2ncvJ4/JommBX1o7V0QiU0soF6BV4nAxc8CdnMegG9UlNTS1uXiEhYiImO4oW+lxETHcUrX6/ndG4e/3p90zIPjlDpnloMNDGzhmYWB/QDZl7oyfQ8DRGJRNFRxrN3tKJ/u/qs3HqQnLyy72wJekvDzKYA3YCaZpYN/NY5N8HMhgNfANHAROfcj6W4hloaIhKRoqKMP9/WklO5ecRGl327wJzzb7d/Wlqay8jI8LoMEZGwYmZLnHNpRb0XKt1TIiISBnwZGmbWy8zGHjwYWo9JFBEJd74MDQ2Ei4gEhi9DQ0REAsOXoaHuKRGRwPBlaKh7SkQkMHwZGiIiEhi+XqdhZruBsxvNJwDn6q8613s1gT0BKK2snO/v5PW5L+T7JflOcT57If/m53tP90Nwz6H7oXRK82/WwDlXq8h3nHMR8QcYW9L3gAyv677Qv5PX576Q75fkO8X57IX8m+t+CNy5S3oO3Q/e/5sV9SeSuqc+vsD3Qlkg6y7tuS/k+yX5TnE+e6H/5rofAnPukp5D90PpBKRuX3dPlZaZZbhzLKWXyKP7QQqK1PshkloaF2Ks1wVISNH9IAVF5P2gloaIiBSbWhoiIlJsCg0RESk2hYaIiBSbQuMCmdmtZjbOzD4ys+u8rke8ZWaNzGyCmc3wuhbxhplVNLM3838u3OV1PYESkaFhZhPNbJeZrSp0vIeZrTWzTDN76nzncM79zTn3ADAIuDOA5UqAldH9sNE5NziwlUqwlfDeuB2Ykf9z4eagFxskERkawCSgR8EDZhYNjAJ6Ai2A/mbWwswuMbNZhf4kFvjqf+R/T8LXJMrufhB/mUQx7w0gGcjK/1huEGsMqhivC/CCc262maUUOtwOyHTObQQws6nALc65vwA3FT6HmRnwV+Az59zSAJcsAVQW94P4U0nuDSCbM8HxAz7+hdy3f7ELkMQ/fkuAMzdA0nk+/yhwDdDbzIYFsjDxRInuBzOrYWajgdZm9nSgixNPneve+AC4w8xeJ3y3HvlFEdnSOAcr4tg5Vz4650YAIwJXjnispPfDXkC/PESGIu8N59xR4L5gFxNsamn8QzZQr8DrZGCbR7WI93Q/yLlE9L2h0PiHxUATM2toZnFAP2CmxzWJd3Q/yLlE9L0RkaFhZlOAdKCpmWWb2WDnXA4wHPgCWANMc8796GWdEhy6H+RcdG/8f9qwUEREii0iWxoiInJhFBoiIlJsCg0RESk2hYaIiBSbQkNERIpNoSEiIsWm0BARkWJTaIiISLEpNEQ8lv98BpGwoF1uRTxgZtM5s712a+Br4BlvKxIpHoWGiDcuAdY457p7XYhISWjvKZEgM7N4YAtwUf7mdyJhQ2MaIsF3MbBQgSHhSKEhEnyXACu8LkLkQig0RIJPoSFhS2MaIiJSbGppiIhIsSk0RESk2BQaIiJSbAoNEREpNoWGiIgUm0JDRESKTaEhIiLFptAQEZFi+x/Nqq8k0zapMAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "rbin = np.logspace(-2,np.log10(2),30)\n", "plt.plot(rbin, hp.nfw(rbin), '-')\n", "\n", "# plt.plot(rbin, hp.esd(rbin)/(1e12), '-')\n", "plt.xscale('log')\n", "plt.yscale('log')\n", "plt.xlabel(r'$r$')\n", "plt.ylabel(r'$\\rho(r)$')\n", "# plt.ylabel(r'$\\Delta \\Sigma (R) [{\\rm h M_\\odot pc^{-2}}]$')" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, '$\\\\Sigma (R)$')" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEMCAYAAADXiYGSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAiO0lEQVR4nO3dd5xV5Z3H8c/vTqXDFDo4MDQH6TMgKIjRICjYK0oTe4kmWTea5mY3MbsxWRNbEKVJLDEYY2+JBSnC0AUUGNrQmcYAAzPDzDz7BxBZMsCdcu+55ft+vXi9mHPvPfPFOfLlOec8zzHnHCIiImfi8zqAiIiEBxWGiIj4RYUhIiJ+UWGIiIhfVBgiIuIXFYaIiPgl1usAgZKSkuLS0tK8jiEiElaWLl2a75xLre61iC2MtLQ0lixZ4nUMEZGwYmZbT/WaTkmJiIhfVBgiIuKXkCwMM+tsZtPMbM4J24ab2RdmNsXMhnuXTkQkOgWtMMxsupntNbPVJ20faWbrzCzHzB4GcM5tcs5NPmkXDjgIJALbg5NaRESOC+YIYyYw8sQNZhYDPAOMAjKAm8ws4xSf/8I5Nwr4EfCLAOYUEZFqBK0wnHNzgcKTNg8Eco6NKMqBV4ErTvH5qmO/LQISqnuPmd1hZkvMbEleXl49JRcREfD+GkY7YNsJX28H2plZsplNAfqZ2SMAZna1mT0HzAaerm5nzrmpzrlM51xmamq1txH75d1Vu8jZexAt/S4i8i2v52FYNducc64AuOukjX8F/hroQEUl5dz78jIAmjeMo1+H5vTv2IL+Z7WgT4fmNE7w+j+ZiIg3vP7bbzvQ4YSv2wM7PcoCQLMGcXz8/WEsyy1i2dZ9LMst4tN1R09v+Qy6t25K/47flkhackPMqus9EZHIYsE87WJmacA7zrlzjn0dC6wHLgJ2ANnAWOfcmjp8jzHAmC5duty+YcOGuocGig8dYfm2Ipbl7mN5bhErcvdxoKwCgKRG8WSe1YKBnZLISksio21T4mK8PtMnIlI7ZrbUOZdZ7WvBKgwzewUYDqQAe4BHnXPTzOxS4PdADDDdOfer+vh+mZmZLlBLg1RWOXL2HmRZbhFLthSRvaWQ3MJDADSMj6Ffx+ZkpSUxMC2Jvh2b0zDe64GciIh/QqIwgi2QhVGdPftLyd5SSPbmQrK3FPH17v04B7E+o2e7ZgxMa8GgTskM7JxE08S4oOUSEakJFYYH9pceYenWomMFUsjKbcWUV1bhMzinXTPO7ZzM4M7JZHVK0oV0EQkZUVUYgbiGUR9Kj1SyPHcfCzcV8OXGApZvK+JIpSPGZ/Rq14zB6cmc2zmZrLQWOoUlIp6JqsI4zusRxpkcLq9kWW4RCzcWsHBTASu37aOiyhHrM/p2aM75XVMY2jWFPu2bE6uL6CISJCqMMHCovIIlW4pYuKmABTn5rNpRjHPQJCGWwenJDO2WytAuKZyl23hFJIBOVxg69xEiGsbHMqxbKsO6HZ2hXlRSzoKNBczLyeOLDfl8tHYPAO1bNGBo11SGdk1hSHoyzRvGexlbRKKIRhhhwDnH1oJDfLHhaHks3FjAgbIKfAZ9OzTnwu4tubBHS3q2barRh4jUSVSdkgrVi971qaKyipXb9zF3fT6frdvLyu3FALRsknCsPFI5r0sKTXT7rojUUFQVxnGRNMI4k7wDZcxdn8cn6/Yyd30eB0oriIsxstKS/jn6SE9tpNGHiJyRCiOKVFRWsXTr0fWvPv1mL+v2HACgY1JDRmS0YkTP1gw4qwUxPpWHiPwrFUYU27HvMJ9+s5dPvtnLvA35lFdWkdQonot6tGREz9YM7ZpCYlyM1zFFJESoMASAg2UVfL4uj4/W7uaTb/ZyoLSCBnExDOuWwoiM1nynR0taNNJdVyLRLKoKIxoueteH8ooqFm0u4KM1e/ho7W727C8jxmcMTEtiVK/WjDynNS2bJHodU0SCLKoK4ziNMPxXVeX4akcxH63dzQerd7MxrwSfwaBOyVzWuw0jz2lNSuNqn4orIhFGhSF+c86xfs9B3l21k3dW7WJT/tHyGJyezOjebbmkZ2uSdNpKJGKpMKRWnHN8s/sA767axTurdrKl4BAxPmNIejKje7fhkp6tNdNcJMKoMKTOnHOs2bmfd7/axburdpFbeIi4GOPC7i25un87LuzRkoRY3W0lEu5UGFKvnDt6zeOtFTt5c+VO8g6U0axBHJf1bsPV/dox4KwWmiQoEqZUGBIwFZVVzMvJ543lO/hwzW5Kj1TRMakhV/Zrx1X92tEppZHXEUWkBqKqMHRbrXcOllXwwerdvLF8Ows2FuAc9OvYnKv7tWNMn7a63iESBqKqMI7TCMNbu4oP8+aKnbyxbAfr9hwgPtbHpee05saBHRnUKUmnrERClApDPHP8Yvmfs7fxtxU7OFBaQaeURtyQ1YFr+rcntYnmd4iEEhWGhITD5ZW899UuXs3OJXtLEbE+4+KzW3HjwA4M7ZqqBRFFQoAKQ0JOzt6D/Dk7l9eX7aCwpJx2zRtwXWZ7rs/sQNvmDbyOJxK1VBgSssorqvh47R5ezc5lXk4+BozIaM34IWcxuHOyrnWIBJme6S0hKz7Wx2W923BZ7zZsKzzEy4tzeXVxLh+s2U33Vk0YP+QsrurXjobxOlRFvKYRhoSc0iOVvLVyJ7MWbGHNzv00TYzl+swOjB+cRsfkhl7HE4loUXVKSvMwIodzjqVbi5i1cCvvf7WLSue4sHtLJgxJY2iXFHy6SC5S76KqMI7TCCOy7NlfykuLcnl5US75B8vonNKIW8/vxLUD2uuJgSL1SIUhEaOsopL3v9rNjPmbWbm9mORG8Uwcksa4wWdpJrlIPVBhSMRxzrF4cyHPzd3EJ9/spUFcDDdkdWDy+Z3okKTrHCK1pbukJOKYGYM6JzOoczLrdh9g6txN/OnLrcz+ciuje7fhjmGd6dm2mdcxRSKKRhgSMXYVH2b6vM28sngbB8sqGNo1hTuHpXNeF83nEPGXTklJVCk+fISXF+Uyff5m8g6UcU67pjxwUTcuPrulikPkDFQYEpXKKip5Y9kO/vj5RrYWHKJn26Z876KujMhopeIQOQUVhkS1isoq3lyxk6c+2cCWgkNktPm2ODSXQ+T/U2GIcLQ43lq5k6c+yWFzfglnt2nKAxd1YURGaxWHyDFRVRia6S1nUlFZxdurdvLUP3LYlF9Cj9ZNeOCirlzSU8UhElWFcZxGGHImlVWOt1fu5MlPNrApr4TurZrwwxHd+K6ucUgUO11h+IIdRiRUxPiMK/u14+PvX8AfbuzLkaoq7pi9lGunLCR7S6HX8URCjgpDol6Mz7iibzs+enAYv766F9sKD3HdlIXcNiub9XsOeB1PJGTolJTISQ6XVzJ9/mamfL6RkrIKrunfnu9/t5ueBChRQdcwRGqhqKScZz/LYdaCrWAwcUga9wxP1yKHEtFUGCJ1sL3oEE98vIG/Lt9O44RY7h6ezqQhnWgQr2XVJfLoordIHbRv0ZDfXd+H9x8YysC0JH7zwTou/O1nvLliB5H6Dy6R6qgwRPzUo3VTpk3M4s93nEtKk3geeHUF101ZyOodxV5HEwkKFYZIDQ3qnMyb957Pf1/di835JYx5eh4Pv76K/INlXkcTCSgVhkgtxPiMGwd25NOHhjP5vE7MWbqdCx//jBe+2ER5RZXX8UQCQoUhUgdNE+P46egMPvz+MAakteCX737NyD/M5dN1e72OJlLvVBgi9SA9tTEzJw1k+sRMnINJM7K5dWY2m/IOeh1NpN6oMETq0Xd6tOLDB4fx40t7sHhzIZf8fi6/+eAbSo9Ueh1NpM4irjDMbIyZTS0u1p0r4o34WB93DEvn038bzuV92vHsZxsZ8cRcPl+f53U0kTrRxD2RAFuwMZ+fvrGaTfklXN6nLT8bnUFqkwSvY4lUSxP3RDw0JD2F9x8cyoMXd+WD1bu56Hef8fKiXKqqIvMfaxK5VBgiQZAQG8ODF3fj/QeHktG2KT9+4yuue24h63ZrNVwJHyoMkSBKT23MK7efy2+v68OmvINc9uQXuiguYUOFIRJkZsa1A9rzjx8O58p+314Un6uL4hLiVBgiHklqFM9vr+vDK7efS6zPGD99MQ+/vooDpUe8jiZSLRWGiMcGpyfz3gNDuXt4Oq8t2cYlT8zliw0abUjoUWGIhIDEuBh+NLIHr989hAbxMYybtpgfv/EVB8sqvI4m8k8qDJEQ0q9jC9793lDuHNaZVxbncskTc5mfk+91LBFAhSESchLjYnjk0rOZc9dg4mN93PzCIn72t9WUaLQhHlNhiISoAWcl8d73hjL5/E78adFWRv5hLl9uKvA6lkQxFYZICGsQH8PPRmfw2p2DiTHjxqlf8uibqzlUrtGGBJ8KQyQMZKUl8f4Dw5h0XhqzFm5l9FPz9GhYCToVhkiYaBAfw6NjevLybYMoKavgqmfn8/zcTVqTSoJGhSESZoZ0SeGDB4ZxYfeW/Oq9r5kwYzF795d6HUuigApDJAy1aBTPc+MG8NhVvcjeUsjIP3zB39fu8TqWRDgVhkiYMjPGDurIO/efT+umidz24hJ+/uZqLWQoAaPCEAlzXVo24Y17h3Db+Z14ceFWLn96Hl/v2u91LIlAKgyRCJAQG8NPR2fw4q0DKSw5whXPzGfG/M1E6hM1xRsqDJEIMqxbKh88OJTzu6Twi7fXcuvMbIpKyr2OJREiJAvDzDqb2TQzm3PS9kZmttTMRnuVTSTUpTROYNqETH5xeU/m5xQw+ql5LM8t8jqWRICgFYaZTTezvWa2+qTtI81snZnlmNnDAM65Tc65ydXs5kfAa8HIKxLOzIwJQ9KYc/dgzOD65xbqFJXUWTBHGDOBkSduMLMY4BlgFJAB3GRmGdV92MwuBtYCundQxE+92zfn3fuHckG3VH7x9lrue3m5HtAktRa0wnDOzQUKT9o8EMg5NqIoB14FrjjFLi4EzgXGArebWUieThMJNc0axjF1XCYPj+rBB2t2c/nT83UXldSK13/ptgO2nfD1dqCdmSWb2RSgn5k9AuCc+4lz7kHgZeB551zVyTszszvMbImZLcnL0xPLRI7z+Yy7Lkj/57IiVz4zn78s2XbmD4qcwOvCsGq2OedcgXPuLudcunPu1ye9ONM59051O3POTXXOZTrnMlNTUwMSWCScDeqczLvfG8qAs1rw0JxV/PuclZroJ37zujC2Ax1O+Lo9sNOjLCJRIbVJArMnD+L+73ThtSXbufKZ+WzOL/E6loQBrwsjG+hqZp3MLB64EXjL40wiES/GZ/xwRHdmTMpi9/5Sxjw1jw/X7PY6loS4YN5W+wqwEOhuZtvNbLJzrgK4D/gQ+Bp4zTm3po7fZ4yZTS0u1rMCRM7kwu4tefd7Q0lPbcSds5fyxMfrtVy6nJJF6n3ZmZmZbsmSJV7HEAkLpUcq+ckbq3l92XZGZLTif2/oS+OEWK9jiQfMbKlzLrO617w+JSUiISAxLobfXtebn4/O4B/f7OWqZ+azRdc15CQqDBEBjs4Ov/X8Tsy+dSB5B8u4/Ol5fL5et6fLtyKuMHQNQ6RuhnRJ4e37zqdt8wZMmrGY5z7fqCVFBIjAwnDOve2cu6NZs2ZeRxEJWx2SGvLXe4Yw6pw2/Pr9b3jg1RUcLtd8jWhX48I4tmJsTCDCiEjoaBgfy9Nj+/HQJd15e9VOrp2ygB37DnsdSzx0xsIwM5+ZjTWzd81sL/ANsMvM1pjZ42bWNfAxRcQLZsa9F3Zh2oRMcgsOcflT8/hyU4HXscQj/owwPgXSgUeA1s65Ds65lsBQ4Evgv83slgBmrBFdwxCpf9/p0Yq/3XcezRvGccsLi3hN61BFpTPOwzCzOOfcaddD9uc9waZ5GCL1b3/pEe59aRlfbMjnnuHp/NuI7vh81S0JJ+GqTvMw/CmCUCsLEQmMpolxTJ+YxU0DO/LsZxu575VlWrwwivhzDWOcmeUdW85jwrFt55rZL81saeAjikgoiYvx8dhV5/DTy87m/dW7uWHql+QdKPM6lgSBP9cwfg5cCvQFOpnZx8BfgHjgwYAlE5GQZWbcNrQzU24ZwPrdB7jymfms233A61gSYP4UxkHnXLZzLh/4BdAH6OWc+3fn3BeBjScioeySnq157c7BHKms4po/LtDM8AjnT2G0PvYkuwuAVsB259y+wMaqPd0lJRJcvdo34837zqNDUkNunZnN7C+3eh1JAsSfu6TuAHoDvY79agJ8DiwHljvnXg50yNrQXVIiwVVSVsH3XlnOP77Zy63ndeInl51NjO6gCjunu0vqjOsXO+emnrSz9nxbIKM4+oxtEYlyjRJimTo+k1++u5bp8zeztaCEJ2/qRyMtkx4x/BlhmDvDm/x5T7BphCHindkLt/DoW2s4p10zpk/MIqVxgteRxE91fR7Gp2Z2v5l1PGmn8Wb2HTObBUyoj6AiEhnGDU5j6rhM1u85wDV/XMDWAj1bIxL4UxgjgUrgFTPbaWZrzWwzsAG4CXjCOTczgBlFJAxdnNGKl28/l/2Hj3D1swtYtX2f15Gkjmr0iFYziwNSgMOheqeUmY0BxnTp0uX2DRs2eB1HJOptzDvIhOmLKSwp59mb+zO8e0uvI8lp1OmUlJn99PjvnXNHnHO7jpeFmf2+vkLWFz0PQyS0pKc25q93DyEtuRG3zVrCnKXbvY4kteTPKamrj//GzM4zs8YnvDas/iOJSKRp2TSRP995Lud2Tubf/rKSZz7N0VP8wpDfD1A6NtKYDiw3s4uObw5IKhGJOE2OLVx4Rd+2PP7hOh59aw2VVSqNcOLPDdKNzOx14BBHlwXJAGaY2UIgMZDhRCSyxMf6eOL6vrRqmsjUuZvIO1DGEzf0JTFOD/EMB/6MMDoD85xz45xzpc65ZUAWUAzoaXsiUiM+n/HjS8/mZ6MzeH/1bsZPW0zxIT0hIRz48zyMOOfcEydtK3fO/cg5pymcIlIrk8/vxFM39WPFtn1c99wC9uwv9TqSnIHf1zBEROrbmD5tmXlrFjuKDnPtlAXkFhzyOpKcRsQVhlarFQkvQ9JTePn2czlQWsG1UxbouRohLOIKQ/MwRMJPnw7Nee3OwZjB9c8tZHlukdeRpBoRVxgiEp66tWrCnLuG0KxBHDe/sIj5OfleR5KTqDBEJGR0SGrInLsG06FFQybNyObDNbu9jiQnUGGISEg5Pis8o21T7nlpGa9rKZGQocIQkZDTvGE8L902iHM7J/HDv6xk5vzNXkcSVBgiEqIaJcQyfWIWl/RsxX+8vZY//H2D1p/ymApDREJWQmwMz4ztzzX92/PE39fzX+98TZXWn/KMZmqLSEiLjfHx+LW9adoglunzN3OovILHruqFz6e1T4NNhSEiIc/nM34+OoPGCbE89UkO5RVV/Oba3sTG6CRJMEVcYZzwxD2vo4hIPTIzfjiiO/ExPn738XrKK6t44oa+xKk0gibi/ktrprdIZLv/oq78+NIevLNqF/e+tIyyikqvI0WNiCsMEYl8dwxL5xeX9+SjtXu4a/ZSSo+oNIJBhSEiYWnCkDQeu6oXn63P4/YXl3C4XKURaCoMEQlbYwd15PFr+zA/J5+JMxZTUlbhdaSIpsIQkbB27YD2PHFDX5ZsLWL89MXsL9XT+wJFhSEiYe+Kvu14Zmw/Vm3fxy0vLGLfoXKvI0UkFYaIRISR57Rhyi0D+GbXAcY+v4iCg2VeR4o4KgwRiRgXnd2K5ydksjHvIDc9/yX5Ko16pcIQkYhyQbdUZkzMIrfwEDdrpFGvVBgiEnGGdElh+oQsthaW6PRUPVJhiEhEGtIlhWkTsthSUMLNL6g06oMKQ0Qi1nldUpg+MYvN+UdLo7BEd0/VhQpDRCLaecdGGpvzSxj7/JcqjTpQYYhIxDu/awovTMj850ijSKVRKxFXGGY2xsymFhcXex1FRELI0K6pPD/+6C23Y1UatRJxhaHlzUXkVIZ1S+WFY6WhkUbNRVxhiIiczrBuR0caOXkHuWWalhGpCRWGiESdC7qlMnXcADbsPTrSUGn4R4UhIlFpePeWR0tjz0GtcusnFYaIRK3h3Vvy7M39WbtzP5NmZOt5GmegwhCRqHZxRiueuqkfK7btY/KsbD257zRUGCIS9Ub1asP/Xt+HRZsLuWP2EsoqVBrVUWGIiHD0IUz/c3VvvtiQz70vLedIZZXXkUKOCkNE5JjrszrwX1f05O9f7+HBV1dQodL4f2K9DiAiEkrGDU6jrKKKX777NfGxPn57XR9ifOZ1rJCgwhAROcltQztTVlHF4x+uIyHWx2NX9cKn0lBhiIhU594Lu1B6pJKnPskhIdbHf1zeE7PoLg0VhojIKfzgu90oq6hi6txNJMTF8MioHlFdGioMEZFTMDMeGdWDsiOVTJ27icRYHz8Y0d3rWJ5RYYiInIaZ8eiYnpRVVPHkJzk0TozljmHpXsfyhApDROQMfD7jV1f14mBZBY+99w2NE+IYO6ij17GCToUhIuKHGJ/xxA19OVReyU/+9hWNEmK4om87r2MFlSbuiYj4KS7Gx7M392dQpyR+8NpK/r52j9eRgkqFISJSA4lxMbwwIYtz2jblnpeXsSAn3+tIQaPCEBGpocYJscycNJBOyY247cUlLM8t8jpSUIRkYZhZZzObZmZzTth2tplNMbM5Zna3l/lERFo0imf25IGkNklg4oxsvt613+tIARe0wjCz6Wa218xWn7R9pJmtM7McM3sYwDm3yTk3+cT3Oee+ds7dBVwPZAYrt4jIqbRsmsifJg+iYXwM46YtZnN+ideRAiqYI4yZwMgTN5hZDPAMMArIAG4ys4xT7cDMLgfmAf8IXEwREf91SGrI7MmDcM5xywuL2LHvsNeRAiZoheGcmwsUnrR5IJBzbERRDrwKXHGafbzlnBsC3Fzd62Z2h5ktMbMleXl59RVdROS0urRszKxbB7K/9AjjXlhE3oEyryMFhNfXMNoB2074ejvQzsySzWwK0M/MHgEws+Fm9qSZPQe8V93OnHNTnXOZzrnM1NTUgIcXETnunHbNmDExi13FpYyfvpjiQ0e8jlTvvJ64V90qXs45VwDcddLGz4DPgpBJRKRWMtOSmDp+AJNnLmHyrGxmTx5Eg/gYr2PVG69HGNuBDid83R7Y6VEWEZE6G9o1ld/f2JeluUXc89LSiHrUq9eFkQ10NbNOZhYP3Ai8VZcdmtkYM5taXFxcLwFFRGrq0l5t+NWVvfh0XR4P/WUlVVXO60j1Ipi31b4CLAS6m9l2M5vsnKsA7gM+BL4GXnPOranL93HOve2cu6NZs2Z1Dy0iUktjB3XkoUu687cVO/nPd9biXPiXRtCuYTjnbjrF9vc4xUVsEZFwds/wdApLypk2bzPJjeK5/6KuXkeqE68veouIRCwz4yeXnk3RoXJ+9/F6mjeKZ9y5Z3kdq9YirjDMbAwwpkuXLl5HERHB5zP+55re7D98hJ+/uZrmDeIY06et17FqxeuL3vVO1zBEJNTExfh4emx/ss5K4gevreDz9eE5sTjiCkNEJBQlxsXwwsRMurRswl2zl7IsDFe4VWGIiARJ08Q4Zt2aRcumCUyakc36PQe8jlQjKgwRkSBq2eToCrcJsT7GTVvEtsJDXkfyW8QVhibuiUio65DUkBcnD+RweSXjpy+m4GB4LFYYcYWhi94iEg56tG7K9IlZ7Nx3mFtnZlNSVuF1pDOKuMIQEQkXmWlJPD22P1/tKObul5ZRXhHa606pMEREPPTdjFY8dlUv5q7P40evrwrpdacibuKeiEi4uXFgR/IPlvHbj9aT2iSBH196tteRqhVxhaGZ3iISju69sAt7D5Qxde4mUhsncPuwzl5H+hcRd0pKF71FJByZGY+O6cllvdrwq/e+5o3l272O9C8iboQhIhKuYnzG/97Qh8KSch76yyqSGiVwQbfQedx0xI0wRETCWUJsDM+NH0DXVk24+09LWbltn9eR/kmFISISYpomxjFrUhbJjeOZNDObTXkHvY4EqDBEREJSy6aJvHjrIAwYP30xe/eXeh1JhSEiEqo6pTRixqQsCkvKmTAjm/2lRzzNE3GFobWkRCSS9G7fnCm3DGDDngPcNXupp7PBI64wdFutiESaYd1S+c21vVmwsYCH5qz0bDa4bqsVEQkDV/dvz+79pfzmg3W0bprIIx7MBldhiIiEibsvSGfXvlKem7uJ1s0SmXRep6B+fxWGiEiYMDP+4/Ke7Nlfyn++s5bWTRMZ1atN0L5/xF3DEBGJZDE+48mb+tG/Ywse+PMKFm8uDNr3VmGIiISZxLgYXhifSfsWDbhtVjYbgvRscBWGiEgYatEonlmTBpIQF8OE6YvZXRz4iX0RVxiahyEi0aJDUkNmTMyi+PARJs5YHPCJfRFXGJqHISLR5Jx2zZgybgA5ew8GfGJfxBWGiEi0Gdo1OBP7dFutiEgEuLp/e3YVl/L4h+to3SyRR0bV/8Q+FYaISIS4Z3g6u4tLWbtzP+UVVcTH1u9JJBWGiEiEOD6xr7LK1XtZgApDRCSixPiMGJ8FZN+66C0iIn5RYYiIiF9UGCIi4peIKwzN9BYRCYyIKwzN9BYRCYyIKwwREQkMFYaIiPjFnPPmYeKBZmZ5wNZjXzYDTnVR41SvpQD5AYhWX073Z/J637X5fE0+4897a/MzP91rOh6Cuw8dD3VTl5/ZWc651Gpfcc5F/C9gak1fA5Z4nbu2fyav912bz9fkM/68tzY/cx0Pgdt3Tfeh48H7n1l1v6LllNTbtXwtlAUyd133XZvP1+Qz/ry3tj9zHQ+B2XdN96HjoW4CkjtiT0nVlZktcc5lep1DQoOOBzlRtB4P0TLCqI2pXgeQkKLjQU4UlceDRhgiIuIXjTBERMQvKgwREfGLCkNERPyiwqgFM7vSzJ43szfNbITXecRbZtbZzKaZ2Ryvs4g3zKyRmc069vfCzV7nCZSoKwwzm25me81s9UnbR5rZOjPLMbOHT7cP59zfnHO3AxOBGwIYVwKsno6HTc65yYFNKsFWw2PjamDOsb8XLg962CCJusIAZgIjT9xgZjHAM8AoIAO4ycwyzKyXmb1z0q+WJ3z0p8c+J+FrJvV3PEhkmYmfxwbQHth27G2VQcwYVFH3TG/n3FwzSztp80Agxzm3CcDMXgWucM79Ghh98j7MzID/Bt53zi0LcGQJoPo4HiQy1eTYALZztDRWEMH/EI/YP1gNtePbfx3A0R9+u9O8/37gYuBaM7srkMHEEzU6Hsws2cymAP3M7JFAhxNPnerY+CtwjZn9kfBdTuSMom6EcQpWzbZTzmh0zj0JPBm4OOKxmh4PBYD+4RAdqj02nHMlwKRghwk2jTCO2g50OOHr9sBOj7KI93Q8yKlE9bGhwjgqG+hqZp3MLB64EXjL40ziHR0PcipRfWxEXWGY2SvAQqC7mW03s8nOuQrgPuBD4GvgNefcGi9zSnDoeJBT0bHxr7T4oIiI+CXqRhgiIlI7KgwREfGLCkNERPyiwhAREb+oMERExC8qDBER8YsKQ0RE/KLCEBERv6gwRILIzO40s91mttLMNprZeK8zifhLM71FgsjMngG+cs5NMbOBwHvOuRSvc4n4QyMMkeDqBaw79vvNQLmHWURqRIUhEly9gHXHntp4H/ATj/OI+E2npESCxMw6cHRUsZqjT2lbBVzs9D+hhAmNMESCpzcw1znXF+gG9AAGe5pIpAZUGCLB0wtYDuCcKwJeBi7zNJFIDagwRILnn4VxzNvApR5lEakxXcMQERG/aIQhIiJ+UWGIiIhfVBgiIuIXFYaIiPhFhSEiIn5RYYiIiF9UGCIi4hcVhoiI+OX/ANwGhwhVTRgJAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "rbin = np.logspace(-2,np.log10(2),30)\n", "plt.plot(rbin, hp.sigma_nfw(rbin), '-')\n", "plt.xscale('log')\n", "plt.yscale('log')\n", "plt.xlabel(r'$R$')\n", "plt.ylabel(r'$\\Sigma (R)$')" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, '$\\\\Delta \\\\Sigma (R)$')" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEMCAYAAAAvaXplAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAe9UlEQVR4nO3dd3jV5d3H8ff3nJxMkpBAwh4B2SBTECcqtaBSce9RBbStLX286uhTtdbWqq21fWxdqFVxFidFUbtQUKEaEFmyh4aVACEJgez7+SMBkTBOQs75nZPzeV1XrhN+Z+STix9+vH/jvs05h4iIyP58XgcQEZHIo3IQEZF6VA4iIlKPykFEROpROYiISD0qBxERqSfO6wBNpXXr1q5r165exxARiSrz58/f5pzLOnB7symHrl27kpub63UMEZGoYmYbDrZdh5VERKQelYOIiNSjchARkXpUDiIiUo/KQURE6lE5iIhIPc3mUtbGWvBVITt2VRDnNwJ+H3E+I87vI+A34nx1j3XbA3Xb4+N8tV9+H2bm9a8gItLkYr4cHv73Kj5YUdDo9wf8RrzfR6CuLOL3e0wM+EkM+EiIq31MDPhJjPOTsO97HwkBP8nxflLi40iK95OS4CcpEEdKQu325Pi4fY/xcRroiUh4xHw5/Op7/SjaU0lltaOquoaqGkdldQ1V1Y6qmpra7Xsfq2ufq6iqoeKAx33b67aVV9VQVllNeWUNO3dXUFZZQ1lVNWWV1ZRV1lBeVfvYEPF+H6mJcXVfgYN8HyAtMY60xADpyQEyU+LJSA6QkRxPelKAOL/KRUSCE/Pl0KVVimc/2zlHeVUNeyqqKa2oqnusZndFFbvLq9ldWc3u8qrabeVV7KqooqRs71clJWVVbNtWSklZFbvKqigprzrsz0tLjCMzJZ6WyfF1jwFat0ggq0UC2WkJZKUmkJ2aSHZaAqkJcTpkJhLDYr4cvGRmdYee/GSkxB/159XUOHZVVFG8p5Kduysp3F3BjtIKdu6urHusYMfuSnburiC/pIwVW0oo2FVORVX9EUxiwFdbFKm1xZGdmkiHlkl0zEiiQ0YSHTOSyUgOqEBEmimVQzPi8xlpiQHSEgN0zAjuPc45isuqKCgpI7+4nPyScvJLyigoqfu+uJwVW0qYs3JbvZFJUsC/X1kk0aFlMh0zkujSKpluWS1okaDdSyRa6V9vjDMz0pMCpCcFOCY79bCvLdpTycbCPeQV7iavcA8bd37z/edf7aRoT+W3Xt82LZHu2Sl0z2rxzVd2Cm3TEjXiEIlwKgcJ2t4S6ds+7aDPl5RVkle4hw3bd7OmYFfdVylvLtj4rVFHcryfblkpHJPVgj7t0ujfIZ1+7dNomXz0h9ZEpGmoHKTJpCYG6NMuQJ923y4P5xwFJeWsKSj9Vml8um4Hby3ctO91HTOS6Nc+jf7t0/cVRnZaYrh/DRFB5SBhYGZkpyWSnZbIyO6tvvVcYWkFSzcVs2RTEUs2FrF0UzHvL9267/ms1AT6t0/j2I4tGZGTyeDOGSTF+8P9K4jEHJWDeCojJZ6TerTmpB6t920rKavky80lLNlYtK80PlhZgHMQ5zMGdExneNdMhudkMqxLJunJAQ9/A5HmyZxzXmdoEsOGDXNaCa75KtpTyYINhXy6fgefrtvBorydVFY7zKBXm1SG52RyXNdMRuRk6lCUSAOY2Xzn3LB621UOEo3KKqtZ+PVOPl23g8/W72D+hkJ2V1QD0LttKqf3zub03tkM7pyB36cro0QOReUgzVpVdQ1LNxUzd+12Zi3PJ3dDIdU1jpbJAU7tmcXpvbM5tWeWrogSOYDKQWJK0Z5K5qwq4D/L8/lwRQHbSyvwGQzpnMFpdaOK3m1Tdb+FxDyVg8Ss6hrHorydzFqez7+X57N0UzEAnTKTOHdgB8YPbn/EGwBFmiuVg0idrcVlzFqezzuLN/Px6m3UOOjbLo3xg9szbmB72qUneR1RJGxUDiIHkV9cxtuLNjN94Ua+yCvCDEbkZDJ+UAfG9m+ny2Sl2VM5iBzBum2lTF+4kekLN7FuWynxfh+jemUxfnAHRvdpo8WWpFlSOYgEyTnH4o1FTF+4iRlfbCK/pJzWLRK4fERnrhjRmTa6j0KaEZWDSCNU1zhmryrg+bkbmLUiH78Z3+3flmtGduW4rhm62kmi3qHKQdNniByG32ec1iub03pls2F7KS/M28DfPvuadxZtpnfbVK45oSvnDmpPcrz+KUnzopGDSAPtqahm+sKNPDd3A19uLiYtMY6Lh3XiqpFdPF12VqQxdFhJpIk558jdUMhzn6znvSVbqHaO0X3aMPmMHvTvkO51PJGg6LCSSBMzM47rWjvh39biMl6ct4FnP1nPOcu2MrpPNpPP6MmAjioJiU4aOYg0oeKySp77eD1PfbSOoj2VnN47m8ln9GBgp5ZeRxM5KB1WEgmjkrJKnvuktiR27q7ktF5ZTB7dk0EqCYkwKgcRD+wqr6otiTlrKdxdyaheWUw+oweDO2d4HU0EUDmIeGpXeRVT567nydm1JXFqzyxuH9u73nrbIuGmchCJAKXlVUydu4EnZq+heE8lV4zows3f6UlGitaZEG8cqhw0WYxIGKUkxPGDUd354GejuHpkV1769CtGPfgBz368jsrqGq/jieyjchDxQMvkeO7+Xj/enXwyAzqkc/eMZZz1f3P4aNU2r6OJACoHEU/1bJPK89cPZ8pVQymvquHKp//LxKm5bNhe6nU0iXEqBxGPmRln9mvLP28+hVvH9OLj1dv4zkOzeeC95ewqr/I6nsQolYNIhEiI8/PDUccw62ejGDewPY99sIbTHvyA6Qs30lwuHJHoEZHlYGbjzexJM5tuZmd6nUcknNqkJfKHiwfy1o9OpH3LJCa/spCJU+eTX1zmdTSJIWErBzP7q5nlm9mSA7aPMbMVZrbazG4HcM695ZybCFwLXBKujCKRZFCnlrzxgxP4xVl9mLOqgO/8cTZvfa5RhIRHOEcOzwJj9t9gZn7gEWAs0Be4zMz67veSO+qeF4lJfp8x8ZRuzJx8Mt2zUvjp3zSKkPAIWzk452YDOw7YPBxY7Zxb65yrAF4BzrVaDwDvOucWHOozzWySmeWaWW5BQUHowot4rHtWC169UaMICR+vzzl0AL7e7895ddt+DIwGLjSzGw/1ZufcFOfcMOfcsKysrNAmFfGYRhESTl6Xw8EW4HXOuYedc0Odczc65x4PeyqRCKZRhISD1+WQB3Ta788dgU0eZRGJGgcbRUx6fj7bd5V7HU2aCa/L4TOgh5nlmFk8cCnwd48ziUSNvaOIO87uw4crCzj74Y/IXX/gqT2RhgvnpawvA3OBXmaWZ2bXO+eqgJuA94EvgWnOuaXhyiTSHPh9xoSTu/HGD04gIeDjkinzmDJ7jQ4zyVHRlN0izUhxWSW3vbaId5dsYXSfNvzhooGkJwe8jiURrNlO2W1m48xsSlFRkddRRDyXlhjg0SuG8MtxfflwZT5n/3kOi/J2eh1LolDUl4NzboZzblJ6errXUUQigpnx/RNzmHbDSJyDCx+by9S563WYSRok6stBRA5ucOcM3v7xSZzUozV3TV/KTS9/TklZpdexJEqoHESasYyUeJ66ehi3jenNe0u28L2/fMyyTcVex5IooHIQaeZ8PuMHo7rz0oQRlJZXcd6jH/Nq7tdHfqPENJWDSIwY0a0VMyefzNAuGdzy2iLuf3c5NTU6DyEHp3IQiSGtWyTw3HXDuXxEZx7/cA0/emkBeyqqvY4lEUjlIBJjAn4f947vzx1n9+G9pVu4dMpc8ks0eZ98m8pBJAaZ1d5V/cSVQ1m5dRfnPfIJy7foRLV8I+rLQTfBiTTemf3a8uqNI6mqqeHCx+Yya0W+15EkQkR9OegmOJGj079DOm/96EQ6ZSZz/bOfMXXueq8jSQSI+nIQkaPXLj2J124cyWm9srlr+lJ+NWMp1bqSKaapHEQEgJSEOKZcPYzrTszhmY/XM3FqLrvKq7yOJR5ROYjIPn6fcde4vvz63H58uLKAix6fq2VIY5TKQUTquWpkV56+Zhgbtpdy0RNzySvc7XUkCTOVg4gc1Khe2bwwYQSFpRVc9Phc1hTs8jqShJHKQUQOaUjnDF6ZNJKKqhoueWKuJu2LISoHETmsvu3TmHbjSAJ+H5dOmcuCrwq9jiRhEPXloJvgREKve1YLpt0wkoyUeK586r98smab15EkxKK+HHQTnEh4dMpM5tUbRtIxI4nvP/MZ/1m+1etIEkJRXw4iEj7ZaYm8MmkkPdukMmnqfN5etMnrSBIiKgcRaZDMlHhenDiCwZ1b8pOXP2faZ1o4qDlSOYhIg6UlBph63QhOPKY1t76+iGc+Xud1JGliKgcRaZSkeD9PXTOM7/Zrw69mLOPRD1Z7HUmakMpBRBotIc7PI5cPYfyg9vzuvRU8NWet15GkicR5HUBEoluc38eDFw2korqG37zzJYkBP1ce38XrWHKUVA4ictTi/D7+dMlgyirnc8dbS0gK+LlgaEevY8lR0GElEWkS8XE+Hr1iCCcd05pbXvuCdxZt9jqSHAWVg4g0mcSAnylXD2Volwwmv/I5/1qmG+WiVdSXg6bPEIksyfFx/PXa4+jXPo0fvriAOasKvI4kjRD15aDpM0QiT2pigOeuG063rBQmTs3lv2u3ex1JGijqy0FEIlPL5HhemDCCDi2TuO7Zz/hcs7lGFZWDiIRM6xYJvDjheFq1SOCav37K0k06/BstVA4iElJt0xN5ccIIWiTEcdXTn7Jqa4nXkSQIKgcRCblOmcm8OPF4/D7jiqf+y9c7tCZ1pFM5iEhY5LRO4cUJIyivquGaZz6lsLTC60hyGCoHEQmbnm1SefLqYeTt2MPEqbmUVVZ7HUkOQeUgImE1PCeThy4ZSO6GQm6etpCaGud1JDkIlYOIhN05x7bnjrP7MHPxFu6d+aXXceQgNPGeiHji+pNyyCvcw9MfraN9yySuPynH60iyH5WDiHjCzLjznL5sLtrDb95ZRvv0RMYOaOd1LKmjw0oi4hm/z/i/SwczuFNLJv9tIbnrd3gdSeo0uBzMLMXM/KEI0xiaeE8kuiUG/Dx1zXF0aJnEhKm5rCnY5XUkIYhyMDOfmV1uZu+YWT6wHNhsZkvN7Pdm1iP0MQ9NE++JRL/MlHie+/5w4nzGtc98SkFJudeRYl4wI4dZQHfg50Bb51wn51w2cDIwD7jfzK4MYUYRiQGdWyXz9DXHsa2kguuf+4zdFVVeR4ppwZTDaOfcr51zi5xzNXs3Oud2OOded85dAPwtdBFFJFYM7NSSv1w+mCUbi7jppc+pqq458pskJI5YDs65yqZ4jYhIMM7o04Zfj+/Pf5bnc/eMpV7HiVnBnHO4yswKzCzPzK6p23a8mf3GzOaHPqKIxJorRnThhlO78cK8r3hh3gav48SkYA4r3QWcBQwCcszsn8CrQDzw05AlE5GYdut3e3N672zu/vtS5mklubALphx2Oec+c85tA34FDAQGOOdudc7NCW08EYlVfp/xp0sH0aVVMj98cYGm+Q6zYMqhrZlNMrNTgTZAnnNuZ2hjiYhAWmKAJ68eRmV1DROn5lJariuYwiWYcvglcCxwD7AMGGBm/6q7x+HykKYTkZjXLasFf7l8CCu3lvCzV7/QLK5hEszVSlOcczc55051zmUCOcBDwDZgbKgDioic2jOLn4/tw7tLtvDn/6z2Ok5MOOLEe2Zmzrl9Ve2cywPygJmHeo2ISFObcHIOX24u5o//WkmvtqmM6d/W60jNWlB3SJvZj82s8/4bzSzezE43s+eAa0ITT0Sklpnx2/MHMLBTS26etpDlW4q9jtSsBVMOY4Bq4GUz22Rmy8xsHbAKuAz4o3Pu2RBmFBEBaifpm3LVUFokxDFxai47tA51yARzzqHMOfeoc+5EoAtwBjDYOdfFOTfRObcw1CFFRPZqk5bIE1cNZWtxOT96cQGVmmIjJBo0ZbdzrtI5t1mXsoqIlwZ3zuC+8wYwd+12fvP2Mq/jNEvBTJ9xx2Ge+1OTpmkErecgEpsuGNqRCSfl8NzcDbzy6Vdex2l2ghk5nL/3GzM70cxa7PfcKU0fqWG0noNI7Lp9bG9O7tGaO6cvYcFXhV7HaVaCPqxUN4L4K/C5mZ2xd3NIUomIBCHO7+Mvlw2hTVoiN724gEKdoG4ywZRDipm9DvSidl6lS4CHzOxxIDGU4UREjiQ9OcCjVwxh264K/mfaQt1B3USCKYduwEfOuavqrlxaABwHFAHHhDSdiEgQju3YkjvH9eWDFQU8+oHuoG4KR7xD2jkXOMi2CjObAbQ4yFtERMLuyhGd+WzdDh7650qGdMnghO6tvY4U1Rp0KauZDTKz35nZeuAPgNaOFpGIYGbcd/4Aclqn8JOXF5JfXOZ1pKgWzKWsPc3sLjNbDjwFbAdGOedGADtCHVBEJFgpCXE8duVQSsuruOllrUF9NIIZOSwHzgYudM4Nc8494JxbX/eczvyISETp2SaVe8/rz6frdvDgP1Z6HSdqBVMOFwDrgX+a2fN1N53VOw8hIhIpzh/SkcuGd+LxD9fw7y+3eh0nKgUzt9KbzrlLqL0y6T3gBiDPzJ4B0kKcT0SkUX45rh9926Vx87QvtMRoIwR9Qto5V+qce9E5dw7QB5gHLA5ZMhGRo5AY8PPYlUOocY4fvbSA8qpqryNFlQZdrbSXc26Hc+4J59xpTR1IRKSpdGmVwu8vHMiivCLufedLr+NElUaVg4hItBjTvy0TTsph6twNzPhik9dxoobKQUSavdvG9mZolwxuf30Rq/N3eR0nKqgcRKTZC/h9/OXywSQE/Nz00gLKKnX+4UhUDiISE9qlJ/HgRceyfEsJD7y33Os4EU/lICIx4/Tebbj2hK488/F6Zi3P9zpORFM5iEhMuX1sb3q1SeWW176goKTc6zgRS+UgIjElMeDn4csGU1JWxc9e/ULrPxxC1JeD1pAWkYbq1TaVO87uw4crC3j2k/Vex4lIUV8OWkNaRBrjyuO7MLpPNve/u5xlm4q9jhNxor4cREQaw8x44IJjSU8OMPmVz9lToctb96dyEJGY1apFAg9dPJBV+bu4d+Yyr+NEFJWDiMS0k3tkMemUbrww7yv+sXSL13EihspBRGLez87sRb/2adz2+iK2anlRQOUgIkJ8nI+HLxtMWWUNN09bqMtbUTmIiADQPasFvxzXl49Xb+fJOWu9juM5lYOISJ1LjuvE2P5t+f37K1iUt9PrOJ5SOYiI1DEz7jt/AFmpCUx+ZSGl5VVeR/KMykFEZD8tk+N56OJBrN9eyv3vxu7srSoHEZEDjOzeiutOzOH5eRuYs6rA6zieUDmIiBzELd/tRbesFG59bRHFZZVexwk7lYOIyEEkBvw8dPEgthaXcc+M2Lt7WuUgInIIgzq15IejjuG1+Xn8c9lWr+OElcpBROQwfnJGD/q0S+PnbyxmR2mF13HCRuUgInIY8XE+Hrp4IEV7KrjzrSU4Fxt3T6scRESOoE+7NH46uifvLN7MjEWbvY4TFioHEZEg3HBKNwZ1asmdby0hPwYm51M5iIgEIc7v4w8XD6Sssprb31jc7A8vqRxERILUPasFt43pzX+W5/Nqbp7XcUJK5SAi0gDXntCVETmZ3PP2MvIKd3sdJ2RUDiIiDeDzGQ9eNBDnHLe+tqjZrv2gchARaaBOmcnceU5fPlmznalz13sdJyRUDiIijXDJcZ0Y1SuL+99bztqCXV7HaXIqBxGRRjAzHrjgWBLi/NzSDA8vqRxERBqpTVoid53Tl/kbCnl+3gav4zQplYOIyFE4f0gHTumZxQPvLW9WVy+pHEREjoKZ8dvz+gPwv282n7mXor4czGycmU0pKiryOoqIxKiOGcncNqY3s1cW8MaCjV7HaRJRXw7OuRnOuUnp6eleRxGRGHbV8V0Y1iWDe95eRkFJuddxjlrUl4OISCTw+Yz7LziWPRXV3P33pV7HOWoqBxGRJnJMdgsmj+7BO4s3896SLV7HOSoqBxGRJjTplG70bZfGndOXULS70us4jaZyEBFpQgG/j99deCw7Siu4d+Yyr+M0mspBRKSJ9e+QzqRTujEtN4+PVm3zOk6jqBxEREJg8hk96NY6hdvfWMTuiiqv4zSYykFEJAQSA37uv+BY8gr38OD7K72O02AqBxGREBmek8lVx3fhmU/WseCrQq/jNIjKQUQkhG4d04t2aYnc9toiyquqvY4TNJWDiEgIpSYGuPe8AazK38Ujs9Z4HSdoKgcRkRA7rXc25w3uwKOzVrNiS4nXcYKichARCYM7z+lLamIc//vm4qhYGEjlICISBpkp8fzvWX2Yv6GQv+V+7XWcI1I5iIiEyYVDOzIiJ5P7Zn4Z8TO3qhxERMLEzLj3vAHsqazm3ncie2oNlYOISBgdk92CH5zanbcWboroqTVUDiIiYfbD046ha6tk7nhrMWWVkXnvg8pBRCTMEgN+fjN+AOu37+bRWau9jnNQKgcREQ+c1KM14we157EP17A6f5fXcepROYiIeOQXZ/clKeDnF28uxrnIuvdB5SAi4pGs1ARuH9uH/67bwesLNnod51tUDiIiHrr0uE4M7ZLBve8sY0dphddx9lE5iIh4yOczfnveAErKqrhv5pdex9lH5SAi4rFebVOZcHI3Xp2fx7y1272OA6gcREQiwuQzetAxI4lfvLk4ItZ9UDmIiESApHg/vx7fnzUFpUz5cK3XcVQOIiKR4rRe2Zw9oB1/nrWa9dtKPc2ichARiSB3jetLgt/HndOXeHrvg8pBRCSCtElL5OYzezJn1TbeX7rFsxwqBxGRCHPV8V3o3TaVe2YsY3dFlScZVA4iIhEmzu/jnnP7s6mojEc8mphP5SAiEoGG52Ry/uAOTJm9lrUF4Z+YT+UgIhKhbj+rN4lxfn7596VhPzmtchARiVDZqYn8z3e8OTmtchARiWBXj6w9Of3rt78M68lplYOISATbe3J64849YT05rXIQEYlwe09OPzl7HevCdOe0ykFEJArcflZvEuJ8YTs5rXIQEYkCe09Oz15ZwPtLt4b856kcRESixDcnp5expyK003qrHEREokQ4T06rHEREosjwnEzOq7tzOpQnp1UOIiJR5ud1J6fvDuHJaZWDiEiUyU5N5Kff6cmHKwv4x7LQnJyOyHIws25m9rSZveZ1FhGRSHTNyG+m9Q7FyemwlYOZ/dXM8s1syQHbx5jZCjNbbWa3Azjn1jrnrg9XNhGRaBPn9/Gr7/WjxjnWb2/6cw/hHDk8C4zZf4OZ+YFHgLFAX+AyM+sbxkwiIlFrRLdWfHDLKPq0S2vyzw5bOTjnZgM7Dtg8HFhdN1KoAF4Bzg32M81skpnlmlluQUFBE6YVEYkOCXH+kHyu1+ccOgBf7/fnPKCDmbUys8eBwWb280O92Tk3xTk3zDk3LCsrK9RZRURiRpzHP98Oss0557YDN4Y7jIiI1PJ65JAHdNrvzx2BTR5lERGROl6Xw2dADzPLMbN44FLg7x5nEhGJeeG8lPVlYC7Qy8zyzOx651wVcBPwPvAlMM05tzRcmURE5ODCds7BOXfZIbbPBGaGK4eIiByZ14eVREQkAnl9tdJRM7NxwDig2MxWAelA0SFefqjnWgPbQpOwSRzud4qEz27oZzTk9cG8tjF/54d7TvtD+D8j2Pcc7f5wpOejcZ842r+zLgfd6pxrVl/AlIY+B+R6nbuxv1MkfHZDP6Mhrw/mtY35O9f+EDn7Q0Pec7T7Q3PcJ0K1PzTHw0ozGvlcJAtl7qb47IZ+RkNeH8xrG/t3rv0hNJ/dmM8I9j1Huz8c6flo3CdCktnqmiemmVmuc26Y1zkkMmh/kAPF4j7RHEcOjTHF6wASUbQ/yIFibp/QyEFEROrRyEFEROpROYiISD0qBxERqUflcARmNt7MnjSz6WZ2ptd5xFta31zMLMXMnqv778IVXucJlWZdDg1Zt/pQnHNvOecmAtcCl4QwroRYE+0PWt+8GWrgvnE+8Frdfxe+F/awYdKsy4EGrFttZgPM7O0DvrL3e+sdde+T6PUsTbc/SPPyLMGvcd+Rb1awrA5jxrCK+rmVDsc5N9vMuh6wed+61QBm9gpwrnPuPuCcAz/DzAy4H3jXObcgxJElhJpif5DmqSH7BrWLlHUEFtKM/we72f5ih3HQdasP8/ofA6OBC81MS5c2Pw3aH4Jd31yahUPtG28AF5jZY0TndBtBadYjh0M46LrVh3qxc+5h4OHQxRGPNXR/0PrmseNQa9yXAt8Pd5hwi8WRg9atlv1pf5BDiel9IxbLQetWy/60P8ihxPS+0azLQetWy/60P8ihaN+oTxPviYhIPc165CAiIo2jchARkXpUDiIiUo/KQURE6lE5iIhIPSoHERGpR+UgIiL1qBxERKQelYNIiJjZDWa2xcy+MLM1Zna115lEgqU7pEVCxMweARY75x43s+HATOdca69ziQRDIweR0BkArKj7fh1Q4WEWkQZROYiEzgBgRd1qgjcBv/A4j0jQdFhJJATMrBO1o4Ul1K4etggY7fQPTqKERg4ioXEsMNs5NwjoCfQGRnqaSKQBVA4ioTEA+BzAOVcIvASc7WkikQZQOYiExr5yqDMDOMujLCINpnMOIiJSj0YOIiJSj8pBRETqUTmIiEg9KgcREalH5SAiIvWoHEREpB6Vg4iI1KNyEBGRev4f6U1dkGV1rzIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "rbin = np.logspace(-2,np.log10(2),30)\n", "plt.plot(rbin, hp.esd(rbin)/(1e12), '-')\n", "plt.xscale('log')\n", "plt.yscale('log')\n", "plt.xlabel(r'$R$')\n", "plt.ylabel(r'$\\Delta \\Sigma (R)$')" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "def model(x, log_M200m, c):\n", " log_M200m = np.log10(log_M200m) + 14 # we are looking in units of 10^14\n", " hp = halo(log_M200m, c);\n", " esd = hp.esd(x)/1e12 # remember we are working with pc not Mpc\n", " return esd\n" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[392.98477752 216.4730308 138.35991084 97.34318514 72.89303872]\n" ] } ], "source": [ "print(model(np.linspace(0.2,2.0,5), 13, 6))" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import curve_fit" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "data = np.loadtxt('/home/idies/workspace/Storage/divyar/IAGRG_2022/DataStore/iagrg_dsigma.dat')\n", "x = data[:,0]\n", "y = data[:,1]\n", "yerr = data[:,2]" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "scrolled": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":2: RuntimeWarning: invalid value encountered in log10\n", " log_M200m = np.log10(log_M200m) + 14 # we are looking in units of 10^14\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1EAAAJSCAYAAAAiZk+wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAABcSAAAXEgFnn9JSAACZOklEQVR4nOzdd3hU1dbH8e9KIJSAdEGligKCFAVUQBSwYUERFbDRFBt2vSre+3rV671gw4qKiAIWkKqCiigCUpWiiKKIShUbEZAekuz3jz1hJr3NZFJ+n+c5z2T22eecNTPJZNbsZs45REREREREJHdioh2AiIiIiIhIcaIkSkREREREJA+URImIiIiIiOSBkigREREREZE8UBIlIiIiIiKSB0qiRERERERE8kBJlIiIiIiISB4oiRIREREREckDJVEiIiIiIiJ5oCRKREREREQkD5REiYiIiIiI5IGSKBERERERkTxQEiUiIiIiIpIHSqJERERERETyQEmUiIiIiIhIHiiJEhERKQLMrJKZPWhmM83sNzNzZjY2l8c2M7MDgWO6RzjUsCjI4xURiTYlUSIiIkVDTeDfwInA8jwe+yJwMOwRRVZBHq+ISFQpiRIRESkafgXqOueOBC7N7UFm1g84GXgiUoFFSL4er4hIUVAm2gGIiIgIOOcOAL/k5RgzqwY8DgwDNkcirkjJz+MVESkq1BIlIiJSfA0HdgGP5fcEZvbfwHik2eELS0SkZFMSJSWOmW0KfCBwZnbQzOJyqN/GzJJDjllZWLGWdmZ2X8jz7vJxfA0zG2hmb5jZGjPbExhcv8XM3jGziyN5fFFgZpUDg/NXm9luM9tpZsvM7K6cfvfzcI3DzOxeM1tsZn+GPEdzA9eumskxJ5rZv83sPTP73swSAn+PCWa2yMz+aWbVwxFfOJlZRTM718z+ZWbTzGxjyO/og3k4T2G8LqcAg4FbA606+XVi4LbYvPfpfb7oisTffuj/iVxsc7M5z1lmNinwd73fzPaZ2c9m9qaZnV7wRx854X5PKehzWlzf48PKOadNW4nZgBqAS7e1yuGYz9LVHx3tx1EaNqApsC/0uc/HOQ6me+32AbvTlX0AVIzE8dHegAbA+pBY9wD7Q+6vBKoV8Bpdgd9CznkQ2J7uOWqTyXHPZ/Lc/p2u7E+gQ7Sfx3Rxd8nkPSR1e7CwXhegfKDu2Cz2xwJfAe+GlA0IHNM9j4859fXtnY/nKw5ol8uteX4fb7q6ep8vwlsk/vYDv6PZbQkh534sk+MNeCmTuPamKxsR7ecvi8cf9vf6MDynxfI9PpybWqKkpGkb8vP6wG3LrCqbWV+gM7AF/wcPsCIyoUkqM4sBxuA/OC0pwKnKAF8ANwGNnXMVnHOVgEaB8wOcC4yK0PF5EvgWMU8tGtmcKxaYATTED9A/yzkXD1QE+uK7eJ0AvFmAa3QC3gdqA58ApwLlnHPVAtdpB/wX2JnJ4V8A/wA64P+5V3DOHQZUxn/Y/xM/O9s7ZlYlvzFGyHZgDn6s0eX4DxS5UhivS8BtQBPg9oKcxMyOwr++kL+WqCOBZbnc3ipIrCH0Pl+0hf1v3zlXJ7sN+F9I9TGZnGIAcH3g5ylAk0BcFYFmwLuBfXdYGHogFIf3+jA8p8X5PT48op3FadMWzg0Yiv8n+TvBb0mGZ1G3IrApUOdhgt+cnBTtx1HSN/wHQAe8ATyY+tzn4zxdc9gf+s1jvXAfn494Ux/rg2E41zUhsWX4pg//4T91/xn5OH9F4KfA8VOAmDD/DpwdEt+VeXwOe+Sybh/gvjzGFZtJ2Ybcvm7hel3IpmUGqIL/4PQycEzIdk/gmIGB+xkeSybn6hE4Zidg+Xgdy+Nb73KztcvP482krt7ni/GW37/9HM65JnC+BVnsnxvYvw4ok8n+siHvdxPCEM+DuX3PyMW5Ivpen9/nNBqvc1Hb1BIlJU1o3/6vAj9n9Q3lfUA9/DfOGwNlScDXkQpOwMwa4VsvEoA7CnIu59zcHKqEfnvWLtzHR1n/wO1c51xmrXkTCX5L3y8f578aOBrfReMG51xKPs6RnaUhP9fNzQFm1hO/rtAUM7sgh7p98N/MDjOzc3IblHMuObd1sxDp1wWgGlAJPx5qXcj2aGD/q4H7tXJxrtT3zC+dc87Mzg6MBdtqfozg12Y2JNB6nIFzbr9zbl4ut3CtBVUq3+fNrKqZ/cP8WMQ/LDg2cV5gXMzh0Y4xl/L8t58dM+sIHBe4+0oW1Y4I3K5yziWl3+mcO0jwd6lSQWMKs8J4T0kjl89pTsL6OhdFSqKkpEnt5pHtP1cza4Bvhk7Gd4dJPe5b59z+yIZY6o0G4oE7nXN/Rvhaoa9lbBSOjwgzqwh0Ctz9MLM6zn8VOCtw9+x8XCb1n/G7zrlt+Tg+J51Dfv4pl8e8h2+9jAOmZpVIhSRQsfgPAIUy61whvS4AfwCXZbKNDOz/b+D+9lycKzUhWW1mbwAfARcCFfCtOC3xrT335zPWSCh17/Nm1gv/QfkxfKteVfz4zSOA0/FrhLWKUnh5lZ+//excE7j9G5icRZ2fA7etzSzD8j5mVhZoE7hbZBZ+LsT3lPRy85zmJNyvc5GjdaKkxDC/XkqjwN2VwDf4bxzrmVkV51zouI0n8d1HnnfOfWNmxW52quLIzAYDZwCfOOfGF8Ilu4T8vDoKx0fKcQS/BPsmm3qp++qYWXXn3F+5ObmZlSPY8jbfzI4G/gmcg2/d2A58DrzknMv0H3s25z0CuADftQrgR3x//xw551LMrD9+kPiV+ESql3Pu/ZBrpE+grgt8yCgMBX5dzOxm/Afk1P/PrczsX4GfP3POfeac24vvYpmGmaV+g77QOTcr/f4spL73XQ4k4p/X6c65fWZWFz+OqTNwr5n9N9zPZW4eb7r6pe59PjCm6y387/17+BbHz51zyYEP2d3wXTmLzIf/9Ar6t5/NeSsBvQN33wr8bWTmRfz41mOACWY21Dn3Y+AcTfFLBRyN/7D/VEFiCrOIvtdnJg/PaWbHRuR1LrKi3Z9Qm7ZwbfgP56n9bxsFyr4J3O8UUq9roGwbUB3/YSt1hp4h0X4chfycDQh5zvKzdcnDtY4CdgSe66NDyh9MPV+YH1tVYGvg3J8V1vH4cThZzXaUOvPf7mzqdMzFNVLHsTiymZUMuCik3vF5eAxNQ457JCTuA2Scme/FXJxvP5n//iwE6ufjtYnFJ0qpMZ0f8twnBcpfIR9jfLK43gZyMb4hHK9LyLUy23K6/oBAvVzNzocf9J167m1Ag0zqtA+pUzccz2dBHi+l7H0eP2FA6t9PpuO+8nCuAVk8z2F/vw+5Zlj/9jM5/7Uh52ybQ93bA+8XqfX3hvxObAdeAKrn4drF/r2+oM9pYb3ORXVTS5SUJKnfMm53zq0P/PwV0ALf1WNRYJabZwL7HnDO/WVmLfFdV6D0zdi0Dz84O78S81B3FH5A/L3OuZ9zqlwQgfEbr+O/ETsA3FKIx1cgONtZVuIDW2Zys95H5ZCfs/uWMHRf5SxrZVQt5Oeh+C4dlwNTnXMHzawevltRX+AGM/veOfdMJudJ9Ru+RaASwcc9F7jHObcpD3EBftySmfXDfzN/Ob5F6jn8GLtY/JigwS7w370QFfh1cc41zO/FnXNjgbF5OOTEkJ8HO+c2ZlInIeTngo4XyyAfj7e0vc8/C5TDd0kdWsBzFeb7faqw/u1n4trA7SrnXLavq3PuaTNbh39/OJzg7wP457gy/n9UbltxSsJ7fWZy/ZyGiPTrXDRFO4vTpi1cG35wpQPmhJTdFSgbGbg/JHD/awIzVxH8di4JqBBybDt8s/4q/IfI3fiBkpnOMoP/8HYvvun6AL5bwL/IfCagXNctCRtwVeA5/jL9YyQCLVHAcwS/CRtU2Mdnc97Ux/pgAc9zRUh8x2RT76yQenlZk6Ujab9N7JNJnRh8tyiHn8o2V7+7+A8vd+E/qKQADxfgeYgFJqSLdQxhaoEKuc6G3LxukX5dwr3hJ11w+BbXTGdfxHfhTP1GPazPaz5jDuv7fFHe8FNHp/6etIx2PGF4PGH72w+cr0XI83NzDnUrAm8H6i4L/A3WwLfGnhUoS30vy3bNsVzGVize6wvynBbW61yUN00sISVJZv3dvwrctgysnJ3aR/c2F5yFK/W475xz+0KOvQ8/uP4LfH/zf+G/+X7DzP6dyfWfw/erXoj/Jz4X+A+ZrzGUl7rFWmDGqKfx32IPdpnMjBTm6z0B3By4e4dz7tXCPL6Q7Ar5uWI29UL37cqyVvbn3+ycezt9Bedn63sycLcmadfuyZJz7g/n3JNAd/w/6v/Laaa9bM6VDMxMV/yeC/wnj4JIvy7hlvreN8NlPftim8Dt11F8XkOF+30+W2ZWw8wamlmFnGuH3UWB26+cc0VpTGa+hPNvPyC1xWQ/Oa+R9Dh+nM8PwGnOuY+dcwnOuW3OuY+B0wL7ahKcoKUoKOz3lLw8p5mKwOtcZCmJkhLBzA7DDxiFzP+5Ho9PUqrjuyTNDamT+s81fbP108BRzrnBzrmXnHNP47+hXwr8MzDAOfX6LYEbgBeccwOcc684564FRgCDzKxdfuqWEI/iv/F7GfjezCqFboR0aQgpz003hwzM7DH8N2AA/wi8ZoV2fCHaGvLzUdnUC923NctaGf0S8vP32dT7LuTnBnk4P865L/BfIgBcl5djU5lZb2Bc4G7q7E+TzKxHfs4XBpF+XcIt9b1vWS7qfBnhWHIUoff5zK5jZnatmX2PH1O1HthpZjPNrE3+H0GetQ7cFmRB8iInTH/7cfgeDuBf6+3Z1K0ccp3nM0uiA2XPB+6eWoSmiy+095S8PKe5EY7XuahTEiUlxQn4ViII+efqnEvAr1JfDbgR/+3K3an7A2Nf2qQ/LnDsQpduGtzAt5rT8AvzNQ3Z1Tdw/afTxfV0yP481zWzAYFVz880s0fM7Bcz221m75vZkYE6d5rZz2a238wWmdlx5JKZ9TGz3wqwdczFZVJn0roR/w1Z+i20n39q2WO5fQwhj+Vx/HTG4PthP1GYxxey7/DdJMB/cMxK6r7fXB5mawrUTU2ksmt9sJCf89NKkXqNY7KtldmFfQKVOgvfS0Az/Axmcfh1pKKRSEX0dQknM6uCn40Msp+trijNaBf29/n0AnUn4HsKPIMfF1ke39L6K7DUzC4r+EPJldTxNgnZ1sqlQnq/z618/+0HXIRvNYKc1zFqQnD2x+ym2l4X8nOjLGsVrsJ8T8nLc5pbBX2dizRNLCElReo/+t34JvlQX+EXejPgCefchpB9TQkOgsztAMojA7ehaxy1AxKcc6FvwjjnNpvZVtIu1JqXuqkexQ8KHh54LHcCk81sNtAL3z2wCr7b4WQza5nLrje5GRibnXy1GIVboAteagvSPc65xwvz+MLmnNtrZovwU093x3dVScPMDD+eBfK3TtJsYCBwnJlZFr9PoQn7+kz25yT1Q3yeup8EPsS+if8fNgq4yTnnApNNgB9HMMXMLnXOFdrUuoX0uoRLakKSSBZTJwdafhoH7ka9JYrCeZ9/AD9leId079GrgcFmtgp43cy+dc6tyfMjyJvUD8/hWqi0KL3f5+tvP0Rqt7Mfgfk51A3tqppdi3nocxPNbraHFPJ7Sl6e09wq6OtctEV7UJY2beHY8AtwOmBBJvt64BcifAyIT7cvdcKD5PT7srhOXWAnsDRd+Wr8TDaZHfMFsC6fdQcE4ltCYIB0oPypQPlPpJ0MI3WgeI5TpxaVjQJOLBF4bVMHwt5V2MdH8Xm7JhBzCnByJvt7hzyuM/Jx/s4hx/fNZH/oxBJbCJmYAN86lO0kBPipqlMCxz+ah7guAw4Gjnsp/XXIOP15jzA81xvI5SDxSL8uYfz9uTMQw8ps6pwWqJMIxEUr1pB4Ivo+DxyGn+Xsqhzi+Bh4oxAe7ziCE38cFu3nP5cxR+xvP+T4+oHX0gFDc1G/AsGpzFeQ9WRPiwJ1/iLk/220t8J4T8nHcxrx17k4bFEPQJu2cGzAmsAf6jN5PG5E4Lg1uahbAb/A6D7Szd6DT2a+yOK4z4At+aw7gExmiMN3+XPAQ+nKuwTKr472a5KH1+DB1H8A2dRJfR4cIWuV4FvoUsvvyMe1C3R8lJ+3MvjZx1KTmDMC5TH4RGNnYN8HeXlO09WZHNi/Hb8mStlAeT3SzorXP91xDfEtA9fjv4m0kH318Ml+6joqCUCdXD7mMwkmUKOy+idOxkSqcx6f22r4bi2p26bAuR5LV14p3K9LIf7+pCYkr2RT5/ZAnS+j/fseiCei7/P4BUIPAOVDyg7HdxWtGFI2EPgzF9dtGPI38mA+Hu/pIcfPx3cpLBPYVwk4Gz9dd5H50qygf/vk7n3pwcD+g8ARuYzr2ZDzfoifDj8msLUCPgrZ/0C0n8d0sRfGe32entOCvs4lZVN3Pin2zCye4PikvPbbT51RLNsuHoEBl9MC9S9zzn2drsp+/DoTmSkf2J+fuqnSr9+yI3Cbfv2F1PLqWZy/xDCz+vjui+C/7brXzO7N5pAnXMg4p4Ien01cfQiuUZMfvZxzi3Oq5JxLMrML8TM7NgQ+MbO9+H+s5QPVvgSuLEAsA/AfIk/DTy19IHCN0HWkHnbOjcvk2Nb4liKARDP7G/9FROiaKeuBS5xzv+Uyns/xf6tfAze4wH/s9FxwHSnwrcd5XRfoSzLv9vMPguPmwLcUDEh37cJ4XcIhN2Odisx4qMJ4n8f/rm9zacfC3g/chl+8d16g7BegppnFuKxnNSww59x8M3sEPzPsacByIMnMduEXA08dH3ZnpGLIp0j87QOHxqwNCNz9wDn3ay4PvRc4Ft8lLnU7ENgX+v94AvDfXMZSIt7rC/CcRux1Li6URElJ0IbgJCm5/uca6EfcJqfjzKwsMAn/rV9/59z0TKptIespno/C9zHOT91UWS1ymVW5ZVFeksSk+zmnvv6Vwnx8Vgpt3IFzboOZtcIPou+FHwx9EPgW/2HgOedcfhbITD3/HjPrCgwCrsYPXq6M/xC5IHD+zD4EbMV3MekCnIwfnF8T//u6Cb/22rvAWy4P000753aZ2ZnAnqwSqJC6qYlUOedcdotUhl2kX5eCMrOK+MH2kP175gmB2y8jG1GutCGC7/MBfwA1zKx8ukQqvaPwyVZOCVTojGlLc6ibKefc/5nZPOAm4BSgFv79/Xv8rIrvOOd25OfcERKRv/0QZxL8giPXkx845/aZ2XnAJfjunW3xSbMDNuO70r/mnHs/D7GUlPf6/DynkX6diwXL4f+QSKlmfuX7icClwHXOudFZ1Psffpa5Ji5kMLKZ1cO/oYxwzt2Vj7oDgNeArs65eSF1u+O7JAx0zo0NKW+D/8Bzhyu603OLiBQ5gYk0fsO/17+RTb2Pgd+dc1dlVSdQ71/4KdcXOuc6hzVYEYk6TXEukoVAE/c4fAJ1e1YJVEDqSui3pyu/PWR/fuqKiEghcM79jR/3NsLMjs2sjpndjJ905X+5OGW3wO394YlQRIoSdecTydrj+D7GS4AEM0v/reNi59zPAM65VWb2MnBToO/+AqADfladcc4vOkde64qISKF6GD+RxJJAS9I7+NnamgK34ru1Xu1ymN7czMrh39dnOecWRDRiEYkKJVEiWUsdt9QhsKU3EPg55P7N+AkgrgUux/cZ/jcwLJNj81JXREQKgXMuxcwux7833wW8GNh1EL8GzynOua9ycZ4D+DEzIlJCaUyUiIiISCbMrAZ+/ajfSvIAeRHJOyVRIiIiIiIieaCJJURERERERPJASZSIiIiIiEgeKIkqocxsoJnNNbM/zWyXma0ws3ytZi0iIiIiIkFKokquM4D3gPPwq85PBF43sz5RjUpEREREpJgrkRNLmFlF4GygB9AeaAjEAj8CU4ERzrndUYqtLXAWcBJwMnAkcMA5Vz4Xx5YHhuKnxK6PX7tiFvCAc25LLo7/ANjnnLsk/49ARERERKR0K6lJ1LXA6MDdb4E1+ClKOwKVge+B051zf0QhtneAi9IV55hEBRKoOfjH8Ct+gdaG+GTsT6CDc+6nHM6xGFjtnLs+X8GLiIiIiEiJXWw3Eb9A3lPOuXWphWZ2BPA+vnvb08AVUYhtCbAKWBbYfsvlcffjE6glwNmpLWlmdifwJPAqcHpWB5tZf6AdfpHXXDOz34CKwOa8HCciIiIiUsTVA/Y65+rk9cAS2RKVHTPrACwGDgCHOecSc6jfCTjCOTclmzpVgNuAR5xzKXmMx5FDS5SZlQX+AKoCJzrnvky3fxXQCmjnnFuRyfEXAW8DNzrnXstjfH+XK1eucuPGjfNymIiIiIhIkfbTTz9x4MCBXc65w/J6bElticrOqsBtOaAGvmtcpgLJy1vAkWaW5Jx7J5M6lYAPgQ7ABmB8mOMFOBWfQP2UPoEKmIJPonoAaZIoM+sLjAVucM6Nzce1Nzdu3Lj5t99+m49DRURERESKphYtWrBmzZp89bYqjbPzHR24PYifmCFLzrmDQF9gP/C2mZ0Xuj8wgcUH+ARqPPB62KP1WgduV2axf2W6egCY2WB8AjUgnwmUiIiIiIikUxqTqNsCt7OccwdyquycWwJcACQB08zsLDg00cMMoDN++vBBLnJ9I+sHbrOagW9LunqpY6VexD/eeWZWJ7DVyuliZlbVzBqaWUOgbEpKnnooioiIiIiUaKUqiQq0JF2Db4X6v9we55ybD/QM3H3XzLoD7wDdgOnA1c655LAGm1alwO3eLPbvSVcP4Fb8tO4v4bsspm7LcnG924H1ge3YhISEPIYrIiIiIlJylZokysyOA94ADPiHc25VDoek4Zz7GLgEP47sQ+AcYCbQxzmXFOZw07PUMHLYf4hzrqFzzjLZGubiek8DjQLbuho1auQnZhERERGREqlUJFFmVhe/KG01/EK7z+TzVB8Bn4fcfzQwbirSdgVu47PYXzFwG5YFhJ1zO5xzG5xzG4CDMTGl4tdERERERCRXSvynYzOrCXyMHy/0GnB3Ps8Ti5844lT8DH/JwHtm1iY8kWZrU+C2bhb766arJyIiIiIiEVKipzg3s8r4rnfNgGnA4PxM/mBmMfjFbPsCC4HuwIX47oEfm1lX59w3YQs8o9SuhydmsT+1/OsIxiAiIhJ1zjlK2xqXIpI9M8Msw+iWiCqxSZSZlQPeBdrhu+Fdnp/JH8y/Ii8B/YClwHnOuT3ABDOLw7duzTGzLs6578L2ANJaBOwEGpvZCZmsFXVp4HZmhK4vIiISNcnJySQkJLBr1y4SExOjHY6IFEFxcXFUrlyZGjVqEBsbG/HrlcjufIGudxOArsACoJdzLr/vus8Bg/GL2HZ3zqWOT8I5Nw64HqiFT6SOLVDgWQjE/nzg7vNmdmhsVGAq81bAQudcbmbeExERKTaSk5PZtGkTCQkJSqBEJEuJiYkkJCSwadMmkpMjOWm2V1Jbom4GLg78vA14IYsmvrudc9uyOomZlQEa47vTne2c25m+jnNudKBF6jGgDrAuu8DM7HwyTq8eZ2ZLQ+7/xzn3fro6jwBnAh2BdWa2AGgAnAwkAAOzu66IiEhxlJCQwP79+4mNjaV27drEx8ejCY9EJFRKSgp79uzh999/Z//+/SQkJHD44YdH9JolNYmqFvLzxVnWggfxSVamnHNJZtYLqOic+yubeiPN7D3n3OZcxFYLn/iEsnRlGRbEdc7tN7OuwFDgCvy6VduBccD/5fLaIiIixcquXb4DSO3atalSpUqUoxGRoigmJubQ+8PWrVvZtWuXkqj8cM49iE+QwnGufcC+XNTLVRLjnBsLjC1ALA8ENhERkRLNOXeoC198fFarfIiIeKnvE4mJiTjnIjrZhNrDRUREpEgKnYVPXfhEJCeh7xORnsWzRLZEiYiIiOTG3sSkPNWvGKePTiKiJEpERERKseYPfJSn+huGnx+hSESkOFHbuIiIiIiISB6oJUqkiFNXExGRyFnz8Dlp7u/cd5AOwz49dH/J0G5UqVC2sMMSkSJOn7ZEijh1NRERiZz0XzwdTHYZ9hfFL6fSzzpWpkwZqlSpwhFHHEHbtm3p0aMHF110EWXKFL3YRUoCdecTERERKab69+9P//79ufzyy+nUqRNJSUmMHz+eSy+9lOOOO44vvvgiLNcZO3YsZsaDDz4YlvOJFHf6ekKkiFNXExERycrYsWMzlP3000/cf//9TJo0ia5du7Jo0SLatGlT6LGJlGRqiRIp4lK7koRu2e0XEZHSrXHjxrz99ttcc8017N27l0GDBkU7JJESR0mUiIiISMDGhD3Z3i9OnnzySeLj4/nyyy9ZuHBhmn3vv/8+gwYN4rjjjuOwww4jPj6e1q1b87///Y8DBw6kqdulSxcGDhwIwEMPPYSZHdpSW8Kcc0yYMIG+ffvSpEkT4uPjqVy5MieddBIvvPACKSkphfKYRQqLvrYWERERASYt38zQqavTlF08cjHDLmlJ73b1ohRV/lWpUoVzzz2XKVOmMHfuXE499dRD+6655hr27NlDixYtaNmyJX///TdffPEF//znP5kzZw6zZ88mNjYWgO7du5OUlMSiRYto3bp1mq6BxxxzDAAHDhzgiiuuoFq1ajRv3pwTTzyRbdu2sWTJEoYMGcIXX3yRaddDkeJKSZSIiIiUeuu37WHotNUku7Sz8yU7x9Bpq2nfsDqNasZHKbr8a9OmDVOmTOG7775LU/7SSy9x1llnER8ffEy7du3iiiuuYObMmbz55pv069cPgPvuu486deqwaNEievbsmenkEmXKlGHq1KlccMEFxMXFHSr/888/Oe+88xg3bhyDBg3itNNOi8wDFSlk6s4nIiIipd7byzaTnOIy3Zec4pi0fHMhRxQeNWvWBGD79u1pynv27JkmgQKoXLkyTz31FADvvvtunq5TpkwZevXqlSaBAqhVqxbDhg3L1zlFijK1RImIiEipt2X73hz27yukSMLLBVrW0q8rBbBu3To++OADfvzxR/bs2UNKSsqh+uvWrcvX9b766itmz57Nxo0b2bt3L845du3aVaBzihRFSqJERESk1KtbrWIO+ysUUiThtW3bNgCqV69+qMw5x913381TTz11KGlKLzXxya3ExEQGDBjAhAkTsqyT13OKFGXqziciIiKlXp/29YiNydhaAxAbY8VyYgnwLUMAzZs3P1T29ttvM2LECI466iimTJnCL7/8QmJiIs65QzPzZZVcZWXEiBFMmDCB448/ng8//JDff//90DnXrl2br3OKFGVqiRIREZFSr1HNeIb1asnQqWknl4g1Y3ivlsVyUomdO3cya9YsALp27XqofPr06QC8+OKLXHDBBWmO+fnnn/N1rdRzpiZS4TinSFGmligRERERoHe7ekwf0jFN2fQhHbmsmLZC3XXXXezZs4f27dvToUOHQ+Wpk0zUq5fxcU2aNCnTc6VOGJGUlJTp/vycU6Q4UxIlIiIiEtCgRny294uDn3/+mT59+jBmzBji4+MZM2ZMmv1NmjQB4OWXX07TxW7BggU8/vjjmZ7zyCOPBDjUNS+91HO+9NJLacqnTJnC+PHj8/dARIowdecTERERKaYGDBgAQEpKCn///Tc//PAD33//Pc45jj32WN566y1atmyZ5phbb72VsWPH8sILLzBv3jxatWrFL7/8wsKFC7nrrrt44oknMlznlFNO4fDDD2fKlCl06dKFo48+mpiYGAYNGkTHjh255557mDVrFvfddx+TJ0+mSZMmrFu3juXLl3P33Xdnek6R4kwtUSIiIiLF1Lhx4xg3bhwTJkxgwYIFxMbG0q9fP6ZOncqaNWto165dhmOaNGnCsmXL6NGjB9u2beO9995j9+7djBo1KsuWqPLly/P+++9z1lln8dVXXzF27FjGjBnDDz/8AMBpp53GwoUL6datGz///DMzZ84kLi6OqVOnMmTIkIg+ByLRoJYoERERkWKmoDPdHXfccbz33nt5One7du2YPXt2luc85ZRTmDNnTp7OKVJcKYkSERGRUmtvYlKO98vGBqc+rxinj04ioiRKRERESrHmD3yU7f4Owz5Nc3/D8PMjGY6IFBMaEyUiIiIiIpIHaokSERGRUmvNw+dEOwQRKYaURImIiEippTFOIpIf6s4nIiIiIiKSB0qiRERERERE8kBJlIiIiIiISB4oiRIREREREckDJVEiIiIiIiJ5oCRKREREREQkD5REiYiIiIiI5IGSKBERERERkTxQEiUiIiIiIpIHSqJERERERETyoEy0AxARkfDam5iUp/oV4/SvQEREJC/0n1NEpIRp/sBHeaq/Yfj5EYpERCR3XnvtNcaPH88333zD/v37adKkCXfeeSdXXnlltEMTyZS684mIiIhIVM2ZM4cLL7yQDz74gC+//JK+ffty9dVX8/bbb0c7NJFMKYkSESlh1jx8TpptydBuafYvGdotzX4RKZ7MLNutS5cuaeo2bNgwarHm5I033uCOO+6gffv2HHPMMfzjH/+ge/fuTJo0Kdvjjj32WBo3bpyrawwYMAAzY968eWGIGJ599llatGhBuXLlMjzfhWXp0qVcdNFF1KxZk/Lly9OkSRP+9a9/sXfv3kI932+//cYdd9xBkyZNqFChAtWrV6dt27bcc8892R6X/vVbunQpl1xyCXXq1KFs2bJUr16dM844gylTpuTr8USSuvOJFDMbE/ZkuN+qbtXoBCNFUvoxTgeTXYb9GgclUnL0798/0/JmzZple9yGDRto1KgRp59+etgSi3DasWMH9erVy3L/999/z48//sjtt99eeEEFTJs2jdtuu41q1apx4YUXEh8fn+PzHW5vvvkm/fv3Jzk5mbZt21K/fn2WL1/Of//7X2bOnMmCBQuoXLlyxM+3ZMkSzjvvPHbs2EHz5s258MIL2bVrF2vWrGHEiBE89thjmV4v/es3efJk+vbtS0pKCu3ataNLly5s3bqVefPm8emnn3LvvfcyfPjwfD1XkaD/oiLFyKTlmxk6dXWasotHLmbYJS3p3S7rfzQiIlJyjR07Nsc63333HWXLlo18MGEybtw4li9fzvPPP59lnffeew+AHj16FFZYh7zzzjsATJkyhW7dumVfOQK2bNnCtddeS3JyMq+++ioDBw4E4MCBA1x99dVMnjyZe+65hxdffDGi59u6dSvnnXceBw4cYNq0aVx88cVp9n/xxRdZXjP09UtKSmLIkCGkpKQwceJE+vTpc6jekiVL6NatG4899hiDBw/OdctjpKk7n0gxsX7bHoZOW02yS9uqkOwcQ6etZv22PVkcKSIipV2zZs2KzIfPnLz77rtcf/31jBo1ihNPPDHLejNmzKBKlSp07ty5EKPztmzZAsDRRx9d6NcGnzjv37+fs84661DCA1CuXDlGjhxJxYoVGTNmDAkJCRE933333ceOHTt47LHHMiRQACeddFKW1wx9/b7//nv+/PNPmjVrliaBAujQoQPnnHMOzjlWrFiRq8dTGJRESZGxNzEpT1tp8/ayzSSnuEz3Jac4Ji3fXMgRiYhIcZF+TNSDDz5Io0aNAJg/f36asVQDBgwo0LVSUlKIj4+natWq9OjR41DCEWrQoEGYGdddd12a8tRWiJdeeinNh/n0tm3bxpIlSzj33HMztLBNnTqVk046iQoVKlC7dm369evH1q1bs415yZIlXHTRRdSqVYty5crRsGFDbrrppgzHPfjgg5gZc+fOBaBRo0aHnrfC7BKZmkxkNg6rVq1aNG/enIMHD/LBBx9E7Hzbt29n0qRJVKlShWuvvTZP8ad//cqVK5er46pXr56n60SSuvNJkaFpmbO3ZXv2gzq3bN9XSJGIiEhx16ZNGy655BKmTp1K7dq16d69+6F9p556aoHOvWPHDi688ELmz5/PzJkzKVeuXJqJASZOnMhrr71G8+bNefrppw+Vjx49mltuuYWxY8fSt2/fbK/xwQcfkJycnKEr3/PPP88tt9xCbGwsp59+OjVr1uSTTz7hlFNOoXXr1pme64033mDAgAGkpKTQsWNH6tWrx8qVK3nxxReZNm0a8+bNOzTeqU2bNvTv359Zs2bx+++/c8kll1CpUiUA6tSpk5+nK1/27PG9T6pVq5bp/tRkY9WqVVx99dUROd+iRYs4cOAAZ555JmXLlmXKlCksXLiQgwcP0qxZM3r37k3t2rUzPV/61+/oo4/m6KOP5vvvv2fSpEn07t37UN0lS5bw0Ucf0ahRI0477bQcH0thURIlUkzUrVYxh/0VCikSEREp7nr27EmbNm2YOnUqzZo1y9W4qtyqXr06EyZM4LfffqNhw4bMnz//0L7169dz/fXXU758eSZOnEjFiv5/24gRI7jnnnsYOXIkXbp04bfffgMgNjaWWrVqZbjGe++9R5kyZTj33HMPlW3YsIG7776bcuXKMWvWrEOtKnv37qVnz57MnDkzw3k2b97Mddddh5nx3nvvccEFFwC+Ne2uu+7i6aefpl+/fofG9vTs2ZOePXvSpUsXfv/9d5544okcZz3s0qVLmucgN+bOnZvtbH+pz8nGjRsz3Z9avmHDhlxdLz/n+/bbbwGoXbs2nTt3ZsmSJWmOGTp0KK+99hqXXXZZhvOlf/1iY2MZO3YsPXr0oE+fPjz++OM0btyYX3/9lYULF3LSSSfx+uuvExcXl6vHUxiUREmRkX6q5Z37DtJh2KeH7i8Z2o0qFYrPoNhw69O+HqMX/Jxpl77YGNPEEiJS+jgH+3dGO4q8K18FzMJ2OsviXNu3b6dq1aphu05e1alTh3bt2rFo0SJ++eUXDj/8cC6//HL+/vtvXnjhBVq2bHmo7rPPPktycjI33HADN9xww6HyBg0aZEgEEhMTmT17Np06dUrTcvLqq69y4MABBg8enCYBqVixIs899xzHHXccLt244ldeeYV9+/Zx9dVXH0qgAGJiYhg+fDiTJk1i2bJlLF26lFNOOSVfz0P37t3zPL18Tq1ap59+Om+99RYTJkzg4YcfTpNcLF26lLVr1wKwa9euXF0vP+fbvn07AOPHj6dcuXKMGTOGCy+8kN27d/Pcc88xYsQIrrrqKpo2bUqrVq0OHZfV69e5c2fmz5/PxRdfzPLly1m+fDkAlStX5swzz+TII4/M1WMpLEqipMjQtMzZa1QznmG9WjJ0atrJJWLNGN6rJY1qxkcxOhGRKNi/Ex5tEO0o8u7ejVChathOl9UU50XhW/sTTjiBRYsWsWrVKubPn8/nn39Oz549ufHGG9PUy22LCfhWml27dnHhhRemKV+4cCFAmq5gqZo2bcoJJ5zAypUr05QvWLAAgCuvvDLDMeXKleOyyy7jmWeeYcGCBflOou677758HZedK6+8kv/+979s2rSJiy66iCeeeIL69euzaNEiBg8eTJkyZUhKSiImJnfTH+TnfMnJyQAkJSUxcuRIBg0aBEDNmjV58skn2bRpE1OmTOGxxx7jjTfeOHRcVq/fhAkTGDhwIKeccgoTJ06kRYsWbN26lSeeeIJHHnmEOXPmMH/+/CIzy2Tp/UQqUgz1blePZnUqc+Hziw6VTR/SUetEiYiUYuHsihduJ5xwAgDPPPMMH3/8MfXq1WPMmDEFOueMGTOAjFObp04CUb9+/UyPq1+/foYkKvWYrFqKUstzmpiisMXHxzNz5kwuuOACZs2axaxZsw7tq1+/PnfeeSePPfZYlmOcwnG+1DWjYmJiMk3kBw0axJQpUzJMuJHZ67du3Tr69+9P7dq1ef/994mP918MH3vssYwaNYpff/2VGTNm8Nprr2WYjCRalESJFDMNasRne19ERKSoSE2iZs+eTWxsLG+++WaBZ1ibMWMGTZs25dhjj01TntpVL6vujdnJ6Zj8nDPV8OHD+f777/N0zH333Zfj4r0tW7bk+++/Z/LkySxfvpykpCRat27NFVdcwSOPPAJAixYtcn3NvJ4vNcGsU6dOprPrpe7/448/0pRn9vpNnDiRgwcP0r1790MJVKjevXszY8YM5s2bpyRKJEcH91GVXewg96tti4hIKVK+iu8aV9yUrxLtCApN8+bNiY2NJTk5mQceeKDAazqtWrWKTZs2cffdd2fYd+SRR/LDDz+wcePGDAkWwKZNmzI9Zu3ataxfv54mTZpk2J86ocIRRxyR75hnzZqV54klBgwYkGMSBVChQgX69etHv3790pR/8sknQOZTlofrfKkJ8vbt23HOZUg0U9eUSp29ELJ+/VKnwT/ssMMyjSu1/K+//srT44kkJVFSZJVZP5evyl/P9yn1+DylGWW+3wfHngaVDo92aCIiUhSYhXVsUWmTOmYqKSlyay++8MILh8bO3HbbbQU+33vvvQeQYTwN+KnZ582bx+TJkznzzDPT7Pvhhx/46quvMhzTuXNn5s6dy5tvvsk556Sd4CoxMZHJkycfqpdfhbl+FPh1v1auXEmLFi3o1KlTxM7XsmVLGjVqxPr16/n8888zjBlLfdyhCyZn9fqlTqSROplEesuWLQOy7nYZDVpsV4qsMpv9uJ9mMZvpX+Zj4t+9Bp44Fp5vDzNuh9VT4O9foxukiIhIMVWzZk3Kli3LTz/9dCjRSS91cdn8LMC7cuXKNJMqpE6JXRAzZsygevXqdOzYMcO+gQMHEhcXx/jx4w9NGAGwb98+brvtNlJSUjIcc80111ChQgUmTJjA+++/f6g8JSWF+++/n19++YX27dvne1KJSPrqq68yJMArV67kiiuuwMx47rnnMhzTr18/mjVrxvTp08NyvnvvvReAW2+9lW3bth0qX7FiBU8++SRAmtkWs3r9LrroIgA+++wzXnzxxTT7li5dylNPPQXApZdemiGGaFFLlBRZMdvWZr5j2w9+W/Gav1/9aGjQCRqe6m+raqpvERGRnMTFxdG9e3dmzJhB69atOfHEE4mLi6NTp04MHDgQ4FDikdcZ0Xbv3k3fvn1JTEzk7LPPZvbs2Xz11VeZJj+59euvv7J8+XKuvPJKYmNjM+w/+uijefTRR7njjjvo2rUrXbp0oWbNmixYsICYmBguuOCCDGtF1a9fn5dffpkBAwbQo0cPOnXqdGix3bVr11K7dm3Gjx+f75gj6fbbb2fNmjW0adOGmjVrsmHDBj7//HNiYmIYNWoUXbt2zXDMpk2bWLt2LTt3ZlwaID/nGzx4MHPmzGHy5Mk0bdqUjh07snv3bhYvXkxiYiKDBw8+lPhk9/qdeOKJ3H333TzxxBPcdNNNjBw5kubNm7N161aWLFlCSkoK1113XYYWxmhSS5QUWXt7T6bbgScYevAa3knuSEqlLPoj//UzfPk6TL8enj4enmoJ02+Ala/7fS7jukoiIiLi10m6+uqrSUhI4K233mLMmDFpxu+sWrUKgKuvvjpP573ppptYt24dt912G/fccw8AX375ZYFinTlzJs65DLPyhbr99tuZNGkSbdq0YeHChcyZM4cuXbqwdOlSatSokekxV111FZ999hkXXHAB3333HVOmTGHfvn3ceOONrFixIldjk6Lhqquuonnz5nz11VdMmTKFjRs30rdvX5YtW8bgwYML5XwxMTFMnDiRkSNH0qBBAz799FOWLVtGu3btGD9+PC+//PKhujm9fo8//jjTpk3j7LPP5rfffmP69OmsWbOG008/nTfffJNRo0bl+TFFkqVfdEwklJl927x58+bhaILPq537DtL6odmH7q964Cyq7N8CGxbBxkX+dmfGQaIZVD4SGnYKtlbVOCasixwWtgzPy7/PLtWLEEvO9DsjxVVKSsqhRT6bNm2a6zVvJDxSUlKoUaMG7dq14+OPP871ca+//jr9+vXjhBNOYMmSJezevZuaNWvSrl27Q2Nb8qNHjx589NFHbNu2LcsJCKToKozXL6/vGS1atGDNmjVrnHO5n8YwQN35pPgw8133qh8NJwa+EduxKZBULfS329dnPG7XVlg92W8AlWpDg47BpKpWs2KdVImIiETCl19+yY4dO3j44Ydzfcy6deu46aabiI+PZ+LEiZQrV45y5cpx1FFHsWrVKjZv3ky9evnrdt+5c2cuuOACJVDFVEl7/ZRESfFWtT60qQ9tLvf3/94KGxfDhoW+tWrbDxmP2f07fDvdbwAVawSSqlN9i9XhLUDfdoqISCnXtm1b8tJjKTExkcsvv5zdu3czduzYNFOGd+3alTfeeIN27drRrVs3RowYkedpw1O7BUrxVNJePyVRUrIcdiS0vNRvALv/CHb927gI/liT8Zi9CfDdDL8BlK8a0lLVCeq0gpiMA1hFREQk6L777mPFihVcccUV9O/fP82+xx9/nN27dzNv3jymTZvGuHHjohSlSHgoiZKSrdLh0OJivwHsSYBNi4NdAH/7Bkj3Ldv+HbD2A78BlDsM6p8S7P53RGuI1XgSERGRUCNGjGDEiBGZ7qtTp06m02qLFFdKoqR0ia8Bx/XwG8C+7bBpabD736+rwKVbR+LA37Butt8AysZD/ZN9UtWgExx1IpQpV7iPQ0RERESiRkmUlG4VqkHTc/0GsP9v2Px5MKna+iWkpFvJ/eAe+OlTvwGUKQ912wfXqarbHsqWL9zHIZKNjQl7MtxvVbdqdIIREREpAZREiYQqfxgce5bfABL3BJKqRX7Cil+WQ3Ji2mOS9sOGBX4DiI2Do9oFp1WvdxLExRfu4xAJmLR8M0Onrk5TdvHIxQy7pCW922lhahERkfxQEiWSnbh4aNzNbwAH98GW5YHJKhbClmU+iQqVnOjHXW1aDDwOMWXgyBMDSdWpvitgucqF/lCk9Fm/bQ9Dp60mOd3sWsnOMXTaato3rE6jmkrwRURE8kpJlEhelK0AjTr7DSDpAPyyMrhO1ebP4eDetMekJMGWL/y28CmwWD85xaGk6hSoULXQH4qUfG8v20xySubTEyenOCYt38y93ZsVclQiIiLFn5IokYIoUw4adPDbaf+A5IOw9atgUrVpKSTuSnuMS4atK/22+DnAoE7L4JiqBh2hYvVoPBopYbZs35vD/n2FFImIiEjJoiRKJJxiy0K99n479Q5IToLfvg6uVbVpMezfme4g5+v89jUsfcEXHd4iOKaqQSeoVKvQH4oUf3WrVcxhf4VCikQkf8zs0M8pKSnEaCF0EclGSkpwhuXQ949IUBIlEkmxZfwU6EedCB1vgZRk+P3b4JiqjYth318Zj/vjW7998bK/X7PpoaTKap9UuI9Biq0+7esxesHPmXbpi40xTSwhRZ6ZERcXR2JiInv27KFKlSrRDklEirA9e/xstHFxcUqiREqUmFg4opXfTrkRUlLgz+9DkqpFsOfPjMdtW+u35a9yGDA3rjafpxzH/JTWcKADVFD3P8moUc14hvVqydCpaSeXiDVjeK+WmlRCioXKlSuTkJDA77//DkB8fLxapEQkjZSUFPbs2XPofaJy5chP4KUkSiSaYmKgdnO/nTQYnINt6/yYqo2LfRfAXVszHNYo5ncaxfxOX+bhnn3BT3TR5Fxo2h2q1o/CA5Giqne7ejSrU5kLn190qGz6kI5aJ0qKjRo1arBnzx7279/P1q0Z3w9FREKVL1+eGjVqRPw6SqJEihIzqNXEb+0G+aRq+/rAOlWBcVU7N6U9JOVgcPHfD/8BtY+HJt39AsJHnugTNSnVGtSIz/a+SFEWGxtL/fr1SUhIYNeuXSQmJuZ8kIiUOnFxcVSuXJkaNWoQGxsb8espiRIpysyg+tF+O/FqAP7+7WcefHYUnWNX0zXmK6ranrTH/P6N3xY8AfGHQ5NzfEJ1dBct+isixVJsbCyHH344hx9+OM45nMt86n4RKZ3MLOJjoNJTEiVSzLgq9ZiWchrTUk4jlmRWDapCpQ0fw9oPIeHHtJX3/AFfvu63MuWh0em+y1+T7nDYkdF5ACIiBRCND0siIukpiRIpxpKJJbleR2hyOpz9CGz7EX74ENbOgk1L/JpUqZL2w7qP/MYdfsHfpuf5hOqI1r7VS0RERERypCRKpCSpeQzUvMVPp773L/jxE99C9eMncODvtHV/XeW3ecOg8pGBbn/nQaPToGz56MQvIiIiUgwoiRIpqSpWh1a9/ZaU6Bf6XTvLt1Rt35C27q6tsOI1v5WtCEd3DXb7q3R4VMIXERERKaqURImUBmXi/MQSR3eB7sPgz7Ww9gP4YRZs/gIIGaR9cC+sfd9vGBzVNpBQnQu1W6jbn4iIiJR6SqJEShszOLyZ3zrfCbv/hHWzfQvVj5/CwdDZ/hz8stxvnz4CVeoHW6gangplykXtYYiIiIhEi5IokdKuUi044Uq/JR2ADQv8OKq1s+DvLWnr7twEX7zst7jKcEw330J17NkQH/mF7URERESKAiVRIhJUphwcc6bfznvCrze19kO/bV2Ztm7iLljzrt8sBuqdHFzkt2YTdfsTERGREktJlIhkzgzqtPTb6ffArt/gh498QvXzPEjaF6zrUvyU6puWwCf/hmqN/Ex/TbtD/Q4QWzZqD0NEREQk3JREiUjuVK4Dbfv7LXEvrP8sMDnFR7D7t7R1t6+HpSP9Vr6Kb9lqeh4ccwZUqBad+EVERETCREmUiORdXEXfytS0O6SkwK9f+Zn+1n4Iv32dtu7+nfDNVL9ZLDTo6Lv8NekONRpHJXwRERGRglASJSIFExMDR53ot673w84tgYRqlm+tSj4QrOuS/cQVGxbAR/f7sVNNuvtWqnonQUxs9B6HiIiISC4piRKR8KpSF9pf67cDu+HnuYFFfmfB3m1p6277wW+Ln4UK1f0sf03PhcbdoPxh0YlfREREJAdKokQkcspVguN6+C0lGX5Z4bv8/TAL/liTtu6+v+DriX6LKevXoUrt9letQXTiFxEREcmEkigRKRwxsb7LXr2T4Mx/w/YNgRaqD2HDIkg5GKybctC3YP08Fz68Bw5vEVjk91w4qq3vQigiIiISJUqiRCQ6qjWEU27w2/6d8OMc30K1bjbs25627h/f+m3BkxBfC5qc4xOqxl0hLj4q4YuIiEjppSRKRKKvfBU4vpffkpNgyxd++vS1syBhXdq6e/6EL9/wW2w5OPp03+WvSXeoclR04hcREZFSRUmUSBG3NzEpx/tlY+3Q/YpxxfzPOraMnwa9QUc4+xHY9qPv8rd2ll/M1yUH6yYf8C1X62bD+3dCnVZ+pr9m5/tFgs2yvo6IiIhIPhXzT1siJV/zBz7Kdn+HYZ+mub9h+PmRDKfw1TwGat4CHW+BvX/5bn9rP/C3B3amrfvb136bPxxqNYNWfaBVbz9joIiIiEiYKIkSkeKjYnVodZnfkg/CxsWBNak+8BNVhPrze5jzEMx52M/017ovHHehpk4XERGRAlMSJVLErXn4nGiHUDTFlvXjoY4+Hc75H/y51nf7+/592LIspKILLvD7/t3Q7Dxo1devRRWrt0ARERHJO32CECniiv0Yp8JgBoc389upd0DCT/D1JPj6bdi+PlgvaR98M9Vv8bXg+EuhdR84oo3GT4mIiEiuabEVESl5ajSGrkPh1i9h0GxoNwjKV01bZ8+f8PmL8HIXGHmynz59x+ZoRCsiIiLFjL7iFpGSywzqn+y37sP9LH6rJsIPH6Vd3HfbWj92as5//PipVn2g+UXFdvxUqZvRUUREpJDpP6eIlA5lysFxPfy29y/4drrv7rf585BKIeOnPrjbT5feOnX8VNmohZ5XpX5GRxERkQhTEiUipU/F6tD+Gr/99bMfP7VqYrrxU/vh22l+q1gTWl7qEyqNnxIRESn1lESJSOlW/Wjoch+cfq+f1W/VRJ847dserLN3G3z+kt9qNvVrT7XqA1XrRS/ubGhGRxERkchSEiUiAr51qd5JfksdP/V1YPxUcmKw3ra18Ol//Nawc5EcP6UxTiIiIpGl/7QiIumViYPjLvDb3r9gzTuw6m3YvDRtvRIwfkpERETyTkmUiEh2Klb3U6S3GwR/rQ+sPzXRj6VKlX781PGX+PWnjjxR46dERERKIK0TJSKSW9UbQZd74ZaVcM0n0P5aqFAtbZ292+CLUTC6G4w8CT57AnZsik68IiIiEhFKokRE8soM6rWH85+Eu36Avm/BcRdCbFzaett+8GOnnm4Jr50PK8fD/p3RiVlERETCRt35REQKokwcNDvfb/u2w7fv+PWnNi1JW2/jQr998A9oei606gvHnKHxU4Us/cLDOdEkHSIikhn9dxARCZcK1aDdQL/9tR5WT4ZVEzIZPzXdbxo/VehyWog4PS1ELCIimVF3PimyNibsyfa+SJFWvRGcfk/I+KnBUKF62jqh46eebw+fPa7xUyIiIsWAOeeiHYMUYWb2bfPmzZt/++23hXrdScs3M3TqapJDfj9jzRh2SUt6tyuaC5yK5CgpEX782C/o+8OstOtPhWrQKbj+VIWqhRpiSZe+O9/OfQfpMOzTQ/eXDO1GlQrBLpbqziciUnK1aNGCNWvWrHHOtcjrsfrvIEXO+m17GDotbQIFkOwcQ6etpn3D6jSqGR+l6EQKINfjpxb5LXX8VOu+cMyZGj8VBumTooPJLsN+JU4iIpITdeeTIuftZZtJTsm8hTQ5xTFp+eZCjkgkAlLHTw2aBbetgq7/hOqN09ZJPuAX+p3QF55s6pOqLStAPQhERESiSkmUFDlbtu/NYf++QopEpJBUaxgYP7UCrp2TxfipBPjiZXilGzzfDuY/Dts3RiVcERGR0k5JlBQ5datVzGF/hUKKRKSQmUHddnD+E3DXWug7wY+Lii2Xtl7CjzD3EXimFbx6LqwYC/t2RCNiERGRUklJlBQ5fdrXIzYm86meY2NME0tI6VAmDpqdB73Hw90/QI9noH7HjPU2LYYZt8ETTWBSP/j+Az+BhYiIiESMkigpchrVjGdYr5bEplszJ9aM4b1aalIJKX0qVIW2A2DQh4HxU/+CGsekrZN8ANa8CxMv9+On3r8btizX+CkREZEIUBIlRVLvdvWYPiTtt+7Th3TkMrVCSWlXrSGc/g+4eTlc+ymcdB1UrJG2zr6/YNloeOWMwPipx2D7hmhEKyIiUiIpiZIiq0GN+Gzvi5RqZlC3LZz3uB8/dflEaN4zi/FT/4VnWsOr3WHV25B0ICohi4iIlBRaDENEpLiLLevXk2p6rp9gYs27fv2pjYvS1tu0xG+z/+WnV283CCrXiUrIIiIixZlaokRESpIKVaFtfxj4Adz2NXT7F9Q4Nm2dPX/A/EfhqeNh6rV+7JSIiIjkmpIoEZGSqloDOO0fcPMyP36qzVVpu/ulHITVk/3YqdHd4OtJmtlPREQkF5REiYiUdKnjp3qOhDvXwBkPQOUj09b5ZQVMGwxPHw9zh8Gu36MTq4iISDGgJEpEpDSJrwmd74Lbv4bLxkL9Dmn37/4d5g+Hp1rA1MGwZUVUwhQRESnKlESJiJRGsWWhxcUwaBZcNx/aXJlJV79J8Eo3GH0GfD1ZXf1EREQClESJiJR2R7aBni/4rn7d/i+Trn7LYdq1vqvfvEdh9x9RCVNERKSoUBIlIiJefE047W7f1e/S16DeKWn37/4d5v0PRjSHadf7cVRSIuxNTMrTJiJS2mmdKBERSSu2LBzfy29bv4TPX4ZvpkByoDtfykH4eqLf6raHk2+A4y6EMnHRjVvyrfkDH+Wp/obh50coEhGR4kEtUSIikrUjT4CLX4Q71vg1pyofkXb/lmUw9Rp4uiXMf0xd/UREpFRQS5SIFFt57VZUMU5veflWqZZfc6rT7fDde/D5KNj8eXD/7t9g7n/hs8fh+EvgpOvgqBOjFq7kzZqHz0lzf+e+g3QY9umh+0uGdqNKhbKFHZaISJGlTxQiUmypC1IUxJb1SdLxl8AvK+GLl+GbqcGufsmJsGqC3+qeBCdfD80v8sdJkZX+C4aDyS7Dfn0JISISpO58IiKSP0edCBe/5Lv6df0XVKqTdv+WL0K6+j0Ou/+MTpwiIiJhpq+VRKTYUhekIqJSLTj9H9DptmBXvy1fBPfv+hXmPgKfPQbHXwonX+fHWomIiBRTSqJEpNhSF6QipkwctLzUb1l29XvLb/VO9l39jrtQXf1ERKTYUXc+EREJv0Nd/b6Frv/M2NVv8+cwZZDv6vfZ47BnW3TiFBERyQclUSWEmQ00s7lm9qeZ7TKzFWZ2ZbTjEpFSrtLhcPo9cPtquGSMX1cq1K5f4dNH/AK+79wEW78q1PA2JuzJ9r6IiEhmlESVHGcA7wHnAScAE4HXzaxPVKMSEYFgV79rP4FrP4VWfSAmpBtf8gH46k14+XQYcw58Mw2SD0Y0pEnLN3PxyMVpyi4euZhJyzdH9LoiIlL8KYkqBGbW1szuM7NpZvaLmTkz25+L48qb2UNm9oOZ7TezrWb2qpnVTV/XOXeVc+4p59wy59yPzrnHgVlA70g8JhGRfKvbFnq97Lv6dbkfKtVOu3/zUpgyEJ5uBZ89EZGufuu37WHotNUku7Tj6JKdY+i01azfphYpERHJmpKowvF/wDDgYuDI3BxgZuWBOcADQCXgXWAzMBBYaWaNc3GaqoAGGohI0VS5NnS5F27/Bnq9Ake1S7t/11b49D+Brn5D4NdVYbv028s2k5ziMt2XnOLUGiUiItlSElU4lgAPAz2AOjnUTXU/0DFwbBPnXB/n3MnAXUAt4NXsDjaz/kA7YFR+gxYRKRRl4qDVZTB4TjZd/d6AUafBq+fCt+9AclKBLrll+94c9u8r0PlFRKRk09y/hcA592jofTPLtr6ZlQVuCdwd4pzbHXKuEYEE6TQza+ucW5HJ8Rfhk6frnXMrCxq/iEihqdsW6r4MZ/0HVrwGy8bAnj+C+zct9tthR0H7a+DEARBfI++XqVYxh/0V8nxOEREpPdQSVTSdiu+K95Nz7stM9k8J3PZIv8PM+gJvAzc4516LWIQiIpFUuTZ0uc+Pm+o1Go5qm3b/37/AnIdhxHHw7hD49es8nb5P+3rExmT+hVZsjNG7Xb38Ri4iIqWAkqiiqXXgNqtWpJXp6gFgZoOBscAA59zYvFzQzL7NbANyM/ZKRCQyysRBq94w+FO4dg607J2xq9+Xb8CozvDaebDm3Vx19WtUM55hvVoSm65nQKwZw3u1pFHN+HA/EhERKUHUna9oqh+43ZLF/i3p6mFmdwKPAUOAeWaWOvYq2Tn3Z0SiFBEpTHXb+e3s/8Dy12D5q2m7+m1c5LfD6sJJ18KJ/aFi9SxP17tdPZrVqcyFzy86VDZ9SEda1a0awQchIiIlgVqiiqZKgdusRj7vSVcP4FYgFngJ+DVkW5abCzrnWmS2AT/lOXoRkUiqXAe6DoU7voGLX4YjT0i7/+8t8MmDga5+N8Nvq7M8VYMa8dneFxERyYySqKIptX9J5vPvBvcf4pxr6JyzTLaGEYtSRCSaypSD1n1g8Fy45hM4/lKICelgkbQfvnwdXjoVXjsf1rxX4Fn9REREIIzd+cysX7jO5ZwbH65zFVO7ArdZfSWaOq3U7iz2i4iUHmZQr73f/n7Ez+q3/FXYE9KTeeNCv1WpF5jVL/uufiIiItkJ55iosWTdcpJXpT2J2hS4rZvF/rrp6omICMBhR0DX+6HzXfDtdFj6Ivz6VXD/zs2+q9+84dCqNzGtr4lWpCIiUoyFe2KJVcC7BTi+J9AqPKEUa6sCtydmsT+1PG9z+oqIlBZlykHrvn7h3i3L4PNRsOYdSAl050vaDyvHU3nleCaUbc6ryd35JCWrt1wREZG0wp1EfeWceyi/B5tZQ5REASwCdgKNzeyETNaKujRwO7NwwxIRKWbMoN5Jfvv7Ed/Nb/mrsHfboSodYtfQIXYN36XUp8zaJGh1EcRoyLCIiGQtnP8l/ibr2eRya1/gPKWacy4ReD5w93kzOzQ2KjCVeStgoXMuVzPviYgIvqtft3/6BXx7vgRHtEmz+7iYTcS/M8BPRLHmXUhJiUqYRcHGhD3Z3hcRKe3ClkQ556o6524u4Dlucs5VC1dMRYWZnW9mS1O3QHFcaJmZnZ/usEeAz4GOwDozeztw7JNAAjCw8B6BiEgJUrY8tLkcrpvH7is/YFZy+7T7//gWJvWDlzrBt++UumRq0vLNXDxycZqyi0cuZtLyzVGKSESk6FF/hcJRCzg5ZAM/TXloWa3QA5xz+4GuwH/wLXw9gYbAOOAE59yPhRC3iEjJZUZy3ZO44eAddD8wnA+ST0q7/481MLm/T6a+mVYqkqn12/YwdNpqkl3aeaKSnWPotNWs36YWKRERUBJVKJxzY7NYwyl0G5vJcfuccw84545xzpVzztVxzg1wzunrQBGRMPre1eemg7eza+B8aH5R2p1/rIEpA+HFjvDNVEhJjk6QheDtZZtJTsl8ot3kFKfWKBGRACVRIiIiASmHt4De4+HGxdC8J2nWNv/zO5gyCF7oAKunlMhkasv27Ic2b9m+r5AiEREp2sKSRJnZZWY22cxeNLPW6fbVNLOfw3EdERGRQlG7BfQe55OpFheTJpnathamXgMvnFLikqm61SrmsL9CIUUiIlK0FTiJMrMrgIlARfyscZ+bWb+QKrFAg4JeR0REpNDVbg6XjYWblkCLXqRNpn4IJlNfTy4RyVSf9vWIjbFM98XGGL3b1SvkiEREiqZwtETdCdzrnDvfOdcJuB4YZWZXh+HcIiIi0Xf4cXDZa3DTUjj+UjIkU9OuhZEnw6q3ITkpamEWVKOa8Qzr1ZJYS5tIxZoxvFdLGtWMz+JIEZHSJRxJVBNgWuod59w4oC/wspldFobzi4iIFA2HN4NLx8CQz6HlZaRJphLWwfTr4IWTYdXEYptM9W5Xj+lDOqYpmz6kI5epFUpE5JBwJFH7gTRrOznn3sWvYzQOuDAM1xARESk6ajWFS16BIV9Ay95gIf9OE36E6dfDyJPgqwnFMplqUCM+2/siIqVdOJKob4BT0xc65yYCtwEvhOEaIiIiRU+tJnDJaJ9MteqTNpn66yd45wYY2R6+eqtYJlMiIpK5cCRR44F2me1wzo0G7gU2heE6IiIiRVPNY6HXyzBkGbTqmy6Z+hneuRGebwdfvqlkSkSkBChwEhVYSDbLSSSccyOcc40Keh0REZEir+Yx0GsU3LwcWl+eNpnavh7evQmebwsrX4fkg9GLU0RECkSL7YqIiIRbjcZw8Us+mWpzJVhscN/2DfDezfBcW1g5XsmUiEgxFPYkysxqmdmr4T6viIhIsVOjMfR8AW5ZDm2uSptM7dgI790Cz50IK8ZBUmL04hQRkTyJREvUYUB/M1MrVzFlZlXNrKGZNQTKpqSkRDskEZHirfrR0HOkT6ZOSJ9MbYIZt/qWqRVjlUyJiBQDSnQkM7cD6wPbsQkJCdGNRkSkpKh+NFw0Em5ZASdcDTFlgvt2boIZt/mWqeWvKZkSESnClERJZp4GGgW2dTVq1IhuNCIiJU31RnDR8z6ZOrFfumRqM8y83SdTy8ZA0oGohSkiIplTEiUZOOd2OOc2OOc2AAdjYvRrIiISEdUawoXPwS0r4cT+GZOp9++EZ0+EZa8omRIRKUL06VhERCTaqjWAC5+FW7+EtgMhpmxw399b4P274NkT4IvRSqZERIoAJVEiIiJFRdX60ONpn0y1G5QumfoFPrgbnmnjk6mD+6MVpYhIqackSkRESo29iUkZtuz2R03VenDBU4Fk6pq0ydSurT6ZerYNfP6ykikRkSgok3OVfLEInVdERCTfmj/wUbb7Owz7NM39DcPPj2Q4OataDy4YAZ3vhIVPBRbnDczat+tX+PAfsHAEdLod2vaHshWiGq6ISGkRiZaov4CHnHNaXEhERCQcqtSF85+EW7+C9oMhNi64b9evMOte381v6YtwcF+0ohQRKTXC3hLlnNsOPBTu84qIiBTUmofPiXYIBVPlKDj/CTj1Dlj0NKwYB8mBiSZ2/waz7vMtVp1uh3YD1TIlIhIhGhMlIiKlRsW4MnnaiqwqR8F5j8NtX8FJ10NsueC+3b/DR0PhmdawZCQk7o1amCIiJVXEkigzK2dm9c2scjZ1KgfqxGVVR0RERLJw2JFw3mNw2yo4+UYoUz64b/fv8NH9Ppla/LySKRGRMIpkS9SdwHqgdTZ1Wgfq3BbBOEREREq2w46Ac4f7ZOqUm9ImU3v+gNn/hGdawaJnIXFP9OIUESkhIplE9QTWO+cWZlUhsG8DcHEE4xARESkdKteB7sPgtq/hlCFQJmRM1J4/4eP/g6dbwaJnlEyJiBRAJJOoxsCaXNT7NlBXREREwqFybej+P98y1eHmtMnU3m3w8QM+mVr4NBzYHbUwRUSKq0gmUfFAbr7m2gscFsE4RERESqfKteGc/8LtX2eeTH3yb9/Nb+FTSqZERPIgkknUZqBdLuq1BX6NYBwiIiKlW6XDA8nUauh4K5StGNy3NwE+eRCebgkLRsCBXVELU0SkuIhkEjUbONrMbsmqgpkNwXfly34JeRERESm4SrXg7P/4MVOdbkubTO37C+Y8hHu6JbbwSSoRnM1vb2JSmk1EpLSL5CIYjwJXAU+b2RnAy8BPgAOOAa4DegB/B+qKiIhIYahUC8562LdKLX6OPQtfJN78or22bzuHLRrGgnKVeCXpPF5L7k6HYZ+mOXzD8POjEbWISJERsZYo59xm4EIgIXA7Az/RxHeBn1P39XTObYhUHCIiIpKF+Jpw1kOceuAZXki6kN0uODV6NdvNP8pOYl65O7ksdh4xpEQvThGRIiaiy7E75z4zsyb4VqczgHqBXZuBT4BXnHPbIxmDiIiIZG/Rw5cBl8HeBA5+8QJllo/GAlOgH247eLzsyzxadwmJZz5CSoNToxusiEgREMkxUQA453Y45x5zzp3jnGse2M5xzj2uBEpEwmljwp5s74tI5irGlfFb1dqUPfsh7PZv4NQ70izaG/P7asq/eREVp/WHhJ+iGK2ISPRFPIkKZWaVzaxyYV5TREqHScs3c/HIxWnKLh65mEnLN0cpIpFirGJ1OPNBuHk5HH9p2n3fz4SRJ8NH/4R9O6IRnYhI1EU8iTKzC8zsQzPbCewAdpjZ34GyHpG+voiUfOu37WHotNUkO5emPNk5hk5bzfptapESyZeq9eDSMXDNx3BUyKolKQdhyfPw3InwxWhI1ox9IlK6RCyJMm8M8C5wDlAZ2Imfja9SoOwdMxtrZhapOESk5Ht72WaSU1ym+5JTnFqjRAqq3klw7SfQ6xU4rG6wfG8CfHA3vNQJfvwkevGJiBSySLZE3QYMxC+keyNQxTlX3TlXDagSKPsVuDpQV0QkX7Zs35vD/n2FFIlICWYGrS6Dm5dB139B2fjgvj+/hzcugTcuhT/XRi9GEZFCEskk6jpgL9DZOTfKOXdoCXTn3C7n3CigM7AvUFeKCDOramYNzawhUDYlRdPaStFWt1rFHPZXKKRIREqBuIpw+j/glhXQ5iogpDPJjx/DCx3g/bthT0LUQhQRibRIJlGNgDnOufVZVQjsmxOoK0XH7cD6wHZsQoL+EUrR1qd9PWJjMu8VHBtj9G5XL9N9IlIAhx0BPUfCdfOgQadguUuGZaPhuRNgyUhISoxaiCIikRLJJOpPIDfvnInAtgjGIXn3ND6xbQSsq1GjRnSjEclBo5rxDOvVkth0wytjzRjeqyWNasZncaSIFNiRbWDA+9D7dajWMFi+fyd8dD+8cAp8/z64zMctiogUR5FMoqYD3cysWlYVzKw60A14J4JxSB4F1vba4JzbAByMiSnUmfBF8qV3u3pMH9IxTdn0IR25TK1QIpFnBs0vhCFfwFn/gXKHBff99RNMvALG9YDfVkcvRhGRMIrkp+N/AT8Dn5pZt/Q7A2UfB+rcH8E4RKSUaFAjPtv7IhJhZcpBp1vh1i+h3TVgIR8zNiyAlzrDuzfDrt+jF6OISBiUieC538V31WsLfGxmfwEbA/vqA6l9xJYC76ab5dw5586IYGwiIiISKfE14YIR0P5amP1P+OnTwA4HX74O306HznfCKUOgbPmohioikh+RTKK6hPxs+KQps8E1HTIpU8dpERGR4q52c7hqGqz72CdT237w5Ym7Yc7DsHwsnPUgtOjluwSKiBQTkUyiNOOeiIhIaWcGTc6Gxl1h+Wsw73+wb7vft3MTTBkEn4+Cc4ZB3bbRjVVEJJcilkQ55zbmXEtERERKhdiycPJ1fsHe+Y/BFy9DSpLft/lzeKUbtOwNZ/4bqtSNbqwiIjnQtGsiIiJSeCpUg+7D4KbPoel5afetngTPtYO5/4PEPdGJT0QkF5REiYiISOGreQxcPgH6vQu1jw+WJ+2D+Y/Cc23hq7cgJSV6MYqIZEFJlIiIiETP0V3g+s+gx7MQXytYvutXeOdGGN0VNi6OWngiIplREiUiIiLRFRMLbfvDLSvh1Dsgtlxw369fwWvnwqR+8Nf6qIUoIhJKSZQUGXsTkzJs2e0XEZESpvxhcOaDcPMyaHFx2n1r3oWRJ8HHD8D+v6MSnohIqkhOcZ4jM+sHXA3UBX4BxjnnXo9mTBI9zR/4KNv9HYZ9mub+huHnRzIcERGJlmoN4LKxcPINMGsobF3py5MTYdEz8OWb0O1fcGI/34olIlLIotYSZWbjgaHAeKAfMA4YambjohWTiIiIFCH1T4Fr58DFL0PlI4Ple7fBzNvhpc7w09yohScipVdUWqLM7CagFXCSc25XoHiZmb0DrDCzG5xzL0UjNomeNQ+fE+0QRESkqImJgdZ94LgLYPFzviXq4F6/749v4fWe0KQ7nP0I1Dw2qqGKSOkRrZao24B7QxIoAAL37wRuj0ZQEl0V48rkaRMRkVIkLh663Ae3rIDWl6fd98MseOEU+PBe2PtXdOITkVKl0JMoM6sAHAPMD9z/xMx+CqkyBzjWzMoWdmwiIiJSxB12JFz8EgyeC/U7BMtTkuDzl+DZE2DpS5B8MHoxikiJF42WqNRp1VLnL10JhC4AERe4TS60iERERKR4OepEGPghXDYOqtYPlu/fAbPuhRc6wNpZ4FzUQhSRkqvQkyjn3EF84nR+4P49zrmrQ6pcCqx0zmmJchEREcmaGbToCUOW+anR4yoH9yWsgwl9/Jip37+NUoAiUlJFa0zUI8DjZtYwtNDMjgaGA/+JRlAiIiJSDJUt7xfpvXUltB0AFvLx5ud58NKpMOM22P1ntCIUkRImKkmUc+5d4Bn8THyPmtm1ZjYc+AJ4zDn3XjTiEhERkWKs0uHQ4xm4fgE0Oj1Y7lJgxVg/Xmrh05B0IFoRikgJEbV1opxzjwFn4adZPwsoC5ztnHs8WjGJiIhICVDneOj3Llw+EWocEyxP3AWf/Buebw/fvqPxUiKSb1GdJ9o5txI/PkpEREQkfMyg6bnQ+AxYPgbmDfeTTgDs2AiT+0P9jtD9f3DkCVENVUSKn6i1RImIiIhEXJk4OOVGuPVLOOl6sNjgvk2L4eUuMP0G+Htr1EIUkeJHSZSIiIiUfBWrw3mPwU1L4dhz0u5bNQGeawvzHoXEvdGJT0SKFSVRIiIiUnrUagJXToKrp8PhzYPlB/fCvP/B8+1g1duQopVWRCRrSqJERESk9Gnczc/id8FTULFmsPzvX2D6dTDmTNj0efTiE5EiTUmUZGBmVc2sYWAdr7Ip+jZORERKotgy0G6QX1+q460QGxfc98sKePVsmDwQdmyKXowiUiSFbXY+M/u5AIc751zjcMUiBXY78O/UOwkJCdGLREREJNLKV4Gz/+MTqo8fgO9Clqv8dhp8/z50GAKd74RylaMXp4gUGeGc4rwh4ADLx7FaqKFoeRoYG/h5do0aNY6NXigiIiKFpHoj6PM6bFgEHw2FX1f58uQDsHAEfPkGdPsXnHAVxMRmfy4RKdEi0Z1vBXAHcDRwRC63IyMQh+STc26Hc26Dc24DcDAmRr0+RUSkFGnYCQbPg54vQqU6wfI9f8CMW2HU6bB+QdTCE5HoC+en477ATKAVMAJYBQwDjgf+cM79nt0WxjhERERECiYmBtpcAbesgNPugTLlg/t+Xw3jLoDpN8Lev6IXo4hETdiSKOfcJOfcRfiWpVuAb4EBwGxgi5k9bmZtwnU9ERERkYgrVwm6/dMnUy17p9236i14vj18PRmcRiaIlCZh76flnPvLOfeCc64jvkvfg8DfwF3ACjP7xszuNbN64b62iIiISERUqQuXjIZr58ARrYPle7fBtGvhzUth+8boxScihSqig10C42r+45w7DjgZeA6oAfwPWBbJa4uIiIiEXd12cO2ncPYjUKZCsPzHT+CFU2Dx85CcFL34RKRQFOaMARuBn4Gt+Bn8NFuBiIiIFGl7E5Mybsmwt92N7LtuEcmNugYrH9wLs/8Jr5wBv34dvaBFJOLCOcV5BmZWEegFXAmcAcQCO4FXgPGRvLaIiIhIQTV/4KMcalzLRTHH8UDZ16lhu3zRr1/By12g481w+n0QVzHCUYpIYQt7a5CZxZjZuWb2JvA7MA7oip+57zKgjnPuOufcwnBfW0RERKRwGe+mnMqZBx5nanLnYLFLhkXPwIsd4Ke50QtPRCIibC1RZnYyvsWpD1ALv4DuIuANYJJzbke4riUiIiJSGNY8fE6a+zv3HaTDsE8P3V8ytBtVKpQN3LsMNn8GM++A7Rt80fYN8HpPaH05nP1fiK9RGGGLSISFsyVqCTAE+AO4H2jonDvNOfeyEigREREpjirGlcmwZbefxt3gxiXQ6Taw2GDFVRNgZHtY9bamQxcpASIxuUNj4AFgrZntzeW2JwJxiIiIiBS+uIpw1sNw3Vw4ok2wfG8CTL8O3rgk2FIlIsVSuJMoA8rnY6uQ2clEREREiq0jWvt1pc75H5QNmVzipznwQgdY/JymQxcppsKWRDnnYgqyhSsOERERkSIjtgx0GAI3LYVjzgyWH9wLs/8Fr3SDrV9FLTwRyR8lLyIiIiKRVq0BXDkFLhkDFWsGy39dBaO7+YQqcW/04hORPFESJSIiIlIYzKDlpXDzMmhzZbDcJfuufS+cAj/OiV58IpJrEV1sFw4tuNsOOAIol1U955wW3xUREZGSr2J16PkCtOoNM26H7et9+Y6N8EYvaNXHj6OKr5ntaUQkeiKaRJnZw8AdQHZLdRt+TSklUSIiIlJ6HN0FbloC8x+FRc/6FimAr9+GdR9D92E+oTKLapgiklHEkigzuwf4F5AEzAR+AHZH6noiIiIixU7ZCnDmg3D8JfDerbB1pS/f9xdMvx5WTYQLnoLqjaIapoikFcmWqMHAPqCzc25lBK8jIiIiUrzVaQnXfgJfvAxz/gMHA0to/jzXT4fedSicMsTP9iciURfJiSXqAfOVQImIiIjkQkwsnHIjDFkKx54dLE/aBx8/AKO7wtYvoxefiBwSySTqtwieW0RERKRkqlofrpgEl74K8bWC5b997adD/+ifkLgnevGJSESTqInASWZWLYLXEBERESl5zPw4qSFfwAlXBctdCix5PjAd+ifRi0+klItkEvUgsAaYZmaNI3gdERERkZKpYnW4aCT0nwHVjw6W79gEb1wCU6+FPduiF59IKRW20Ylm9mkmxTFAJ+B7M9sAbMFPZ56ec86dEa5YREREREqURqfBjYvhs8dh0TOQkuTLV0/2LVLn/A9aX67p0EUKSTineOmSzb5YoHFgy0xmiZUUgJmdBtwFtAHqAw855x6MZkwiIiJSAGUrwBkPQIteMONW+GWFL9+3Hd650U+H3uPptC1WIhIR4UyitIBB0VIJ353yLeDp6IYiIiIiYVPneLjmY/hiNMx5ODgd+vr5fjr0LvdBh5shtmx04xQpwcKWRDnnNobrXOFkZnWAe4Hz8dOu7wPWA3Occ/dEIZ62wFnAScDJwJHAAedc+VwcWx4YClyOb136C5gFPOCc2xJa1zn3AfBB4LhHw/kYREREJMpiYuGUG6DZ+fDB3fDDLF+etB8+eRBWT4ULn4WjToxqmCIlVSQnlog6M+sAfAfcDhwE3gOWAjWAO6MU1v8Bw4CL8QlUrgQSqDnAA/hWpneBzcBAYKUm7xARESmFqtaDyyfCpa9B/OHB8t9XwytnwKz74cDu6MUnUkKFc2KJw4BE59z+ApyjPBDnnPs7DPEciW+JKQf0cs5NT7f/pIJeI5+WAKuAZYEtt+tp3Q90DBx/tnNuN4CZ3Qk8CbwKnB72aEVERKRoM4Pje0Hjrn5R3pXjfblLgaUj4bsZcMEIOPassF1yb2JSnupXjAvnCBKR6Avnb/R2YCxwTQHO8SJwNeGJazhQFbglfQIF4Jz7IjcnMbNOwBHOuSnZ1KkC3AY84pxLye58zrk0XessF7PomFlZ4JbA3SGpCVTgfCPMrD9wmpm1dc6tyPGEIiIiUvJUqAYXPgcte8PM2yHhR1++cxO8eSkcfyl0Hw6VamV7mtxo/sBHeaq/Yfj5Bb6mSFESzu58FtjCcZ6CncAv8Nsb2Am8UoDzlMVPzDDBzHpmUacS8CHwEHBVZnXC4FR8QviTc+7LTPanJng9InR9ERERKS4adYYbFsFp/4CYkO+lv5kCz7eDL98Ap4mRRQoi3G2rp5rZqwU5PkxxdMJ34/sEOGhmlwbOXRb4HpjknPs9p5M45w6aWV9gNvC2mV0cmLABADOriO8y2AEYD7wepvjTax24XZnF/pXp6hWImVXFJ20AZVNSsm1cExERkaKmbHno9q/gdOhblvny/Tvg3SHw9dtwwdNQI39Dqtc8fE6a+zv3HaTDsOCSoUuGdqNKBc0OKCVXuJOoYwJbQYTjq5EWgdvfgQX4JCfUMDMb6JybnGMwzi0xswvwydI0M+vhnPs4MH5rBtAZmAgMci5iX+vUD9xuyWL/lnT1UlvIUl+LOKCOmbXBj1tbk8P1bgf+nXonISEhj+GKiIhIkVC7OQz6CJaNgTkPQWJgRMD6z+DFjnD6vdDxljxPh55+jNPBZJdhv8ZBSUkWzt/urmE8V0FVC9z2Aw7gx2m9h5/V7hb8zHxvmNla59zXOZ3MOTc/0J1vBvCumfXCJxrdgOnA1c655HA/iBCVArd7s9i/J109gHbA3JD71we2jUDDHK73NH58G8DsGjVqHJvLOEVERKSoiYmFk6+DZufBB/+AtYFONUn7fWL1zVTo8SzUbRvdOEWKkXCuEzU/XOcKg9jAbRn8RAypXQy3AXeZWX3gUuAecjmOKdD6dAk+afowUDwT6OOcy9sUNXmXOk4sq5auDOPInHPzMivPDefcDmAHgJkdjIkp0TPhi4iIlA5V6kLft2DNu/DhPbA7MLLh929gzJlw0vW+C2C5StmfR0RK7DpRuwK3KcC4TPanJlVd8njej4DPQ+4/6pw7mMdz5Efq44nPYn/FwK0WghAREZGsmUGLnjDkC2g7IFjuUuDzF+GFU+CH2dGKTqTYCFsSZWavmtmgfBx3r5l9mnPNPNkQuP3NOXcgm/2HZ7IvU2YWi5844lT8Ok/JwHuBcUaRtilwWzeL/XXT1RMREZEI2JiwJ9v7xUaFqtDjGRjwAYT22t+5Gd66DCYPhN1/RC08kaIunC1RA8jf7HrNCP8isanTgFezzBdiqhG4zVXLjZnF4Fuv+gIL8bP/XQ1UAT42s+MLFm6OVgVuT8xif2p5juO7REREJH8mLd/MxSMXpym7eORiJi3fHKWIwqBhJ7hhoZ9gIiZkcolvp8Hz7WHl65oOXSQThdKdz8yqF8Z1UjnnVgPrgQrAyZlU6RK4zWrK8EMCSdhL+EkqlgLnOef2OOcmAIPwCdkcMzsuDKFnZRF+zavGZnZCJvsvDdzOjGAMIiIipdb6bXsYOm01yekSimTnGDptNeu3FdMWKfDToXe9H25YAHVPCpbv3wHv3QzjekDCT1ELT6QoingSZWZ3A7+Y2f8FWnQKy6OB22fNrGZIPG2BuwJ3X8rFeZ4DBgMrgO7OudTxSTjnxuFnvKuFT6QiMoudcy4ReD5w93kzOzQ2yszuBFoBC51zyyJxfRERkdLu7WWbSU7JvEUmOcUV79aoVIcf56dDP+8JiKscLN+wAF7oAJ89AcmFMRRcpOiLWFJjZrFmNhqfzJQDHgQ+M7OGkbpmOqOByUB7YK2ZzTCzucBi/EKyo51zU7I7gZmVARrju9Od7Zzbmb6Oc240ftr0KkCdnIIys/PNbGnqFiiOCy0zs/MzOfQR/KQWHYF1ZvZ24PgngQRgYE7XFhERkfzZsj2rVUZS9+8rpEgiLCYGThoMQz6HpiEfR5IPwKf/gVGnw5bl0YtPpIiISBIV6L73MX59ph+BC4HV+ARglZn1j8R1QznnUvBjmIbg10bqhk+olgP9nHPX5eIcSUAv4Azn3F/Z1BsJNHPOLchFaLXwXQxTN/BTkYeW1crkGvvxa3H9B79eVE/8ek/jgBOccz/m4toiIiKSD3WrVcxhf4VCiqSQVDkKLn8Ler8OlUK+I/7jW3jlTPjgHjiwK+vjRUq4sCdRZtYEP3aoC/ApcLJzbiY+gRmBXxD2VTObbGbVsjxRGDjnUpxzLzjnTnTOxTvnKjnnOjnnXs/DOfY55xJyUS9X7fjOubHOOcthG5tNLA84545xzpVzztVxzg3I7bVFREQkf/q0r0dsTObLL8bGGL3b1SvkiApJ8wt9q1Tb0A4vDr4YBSNPgbWzohaaSDSFO4k6AVgCHAOMAs4JLNyKc+6gc+5u4GzgV3wLz2qgdZhjEBEREQmrRjXjGdarJbHpJv2NNWN4r5Y0qpnVUo4lQIWq0ONpGDgLajYJlv+9BSb0gckDsNSFe0VKiXAnUa2Aw4DbnXM3OueS01dwzs0BWgLTgSNREiUiIiLFQO929Zg+pGOasulDOnJZSW2FSq9BBz8depeh6aZDn07lVzrSJ3YuoOnQpXQIdxK1C7jAOfdsdpWcc9udc5cC1+LH94iIiIgUeQ1qxGd7v8QrUw663Ac3LoJ6pxwqtgM7ebTsaCbGPUJ9U6uUlHxhS6KcczHOuarOuY/ycMyrzrnKzrnYcMUhIqXH3sSkDFt2+0VEJExqNYWBH8L5I6DcYYeKT4n5jg/j7qPsqje0SK+UaGWiHYCISH41fyD772w6DPs0zf0NwzNbPUBERPIlJgbaXwNNz+XgzH9Q9oeZAMTbAZh1O2z4BHo8A/E1sz+PSDFUmIvfZmBmzc3sf2a2MZpxiIiIiEg+HXYkey8ey42Jt7HdVQqWfz/TL9K77uPoxSYSIYXeEmVmtYErgKuANoV9fREpOdY8fE60QxARkYAPU05mxYEmPFH2JU6LXe0L9/wBb14K7QfDWQ9DXPbrbZUEee0+XjFOHcOKo0J51cysInAxcDVwBr4FzICd+GlcqhRGHCJSsugfj4hI0fIH1eh/8F7WnL2JCvMeguQDfsey0bB+PvQaDUe2iWqMkZZTV/P01NW8eIpYdz7zzjaz8cDvwHj8GlFJ+OnNLwVqA6siFYOIiIiIFC5HDIntroPr50PtlsEd236AV86ABU9CSoZVcESKlbB/jWtmbfAtTpfjkyQDUoC5wFvAFOfc3yH1wx2CiIiIiETb4cfB4Dkw97+w6FnAQUoSzHnYj5O6eBRUaxDtKMMufVfznfsOppnoaMnQblSpUDb9YVLMhC2JMrN78cnTcfjECWAlPnGa4Jz7NVzXEhEREZFioEw5PxbqmLNg+g3w9xZfvmkJvNgJznscWveFEvSlevqu5geTXYb96o5e/IWzO98wfAL1G/Af4DjnXDvn3AglUCIiIiIl18aEPdnep1Fnv0Bvy8uCZYm74J0bYPIA2PtX5IMUCaNwj4kyfBe+04HOZlY1zOcXERERkSJk0vLNXDxycZqyi0cuZtLyzWkrVqgKl7wCl4yBciFziq15B17sCD/NjXisIuESziTqFOAF4C/gNGAU8JuZTTezS82sXBivJSIiIiJRtn7bHoZOW02yS9tlLdk5hk5bzfptezIe1PJS3yrVsHOwbNev8HpPmDUUDu6PbNAiYRC2JMo594Vz7mbgSKAnMA0/ocRFwNvAH2b2mpmdZZpNQkRERKTYe3vZZpJTXKb7klNcxtaoVFXrQb/34OxHIDYuWL70BXi5C/y2OvzBioRR2Kc4d84lOefec85dBtQBrgcWAZWA/sAsYKuZPWVm7cN9fREREREpHFu2781h/76sd8bEQMdbYPCnUOu4YPmf38Hobn5Gv5SUMEUqEl4RWycKwDn3t3NutHPuNOBo4AHgB/y4qVuBpfiufyIiIiJSzNStVjGH/RVyPkmdlnDdPDjlpmBZciJ8/H8w/kLYuaVgQYpEQESTqFDOuY3OuUecc8cBJ+PHT20jOB26iIiIiBQjfdrXIzYm849ysTFG73b1cneisuWh+zC4+h2ofESwfMMCeKEjrJ5S8GBFwqjQkqhQzrllzrlb8OOnLgL0lyEiIiJSzDSqGc+wXi2JTTfcPdaM4b1a0qhmfN5O2Lgr3LgYml8ULDuwE6ZeA1OvhX07Ch60SBhEJYlK5ZxLds7NcM71iWYcIiIiIpI/vdvVY/qQjmnKpg/pyGW5bYVKr2J1uGwc9HwJ4ioHy1dP9gv0rl9QgGhFwiOqSZSIiIiIFH8NasRnez/PzKDN5X4q9PodguV/b4FxPWD2/0HSgYJdQ6QAlESJiIiISNFUrQEMeB/OeABiygQKHSx+FkafAX98F9XwpPRSEiUiIiIiRVdMLHS+C679BGo2CZb/vhpGnQ5LX9RU6FLolESJiIiISNF35Alw3XxoPzhYlnwAZt0Hb/SCv3+NXmxS6iiJEhEREZHiIa4inP8EXDEZ4g8Plv88F17sAN++E7XQpHRREiUiIiIixUuTs+GmJdDsgmDZvu0wuT9MvxH2/x292KRUUBIlIiIiIsVPfE3o8wZc+ByUDZkNcNVb8FIn2LgkerFJiackSkRERESKJzM4sR/csADqtg+W79gEY8+DOQ9DUmL04pMSS0mUZGBmVc2soZk1BMqmaMYbERERKcpqNIaBs6DL/WCxvsylwIInYcxZ8OcP0Y1PShwlUZKZ24H1ge3YhISE6EYjIiIikpPYMtDlXrhmNlQ/Olj+61cw6jRY9go4V+hhbUzYk+19KZ6URElmngYaBbZ1NWrUiG40IiIiIrlVtx1cvwDaDgiWJe2D9++Ct3rDrt8LLZRJyzdz8cjFacouHrmYScs3F1oMEhlKoiQD59wO59wG59wG4GBMjH5NREREpBgpVwl6PAN9J0DFkC+D1832U6F//37EQ1i/bQ9Dp60mOV3rV7JzDJ22mvXb1CJVnOnTsYiIiIiUTM3OgxuXwLFnB8v2JsDEK+C9W+DA7ohd+u1lm0lOybz7YHKKU2tUMackSkRERERKrsq14YpJcP4IKFMhWL5yPIzqDJuXReSyW7bvzWH/vohcVwqHkigRERGRLOxNTMqwZbdfiigzaH+Nnwr9iDbB8r9+hlfPgbnDIDm8r1/dahVz2F8h2/1StJWJdgAiIiIiRVXzBz7Kdn+HYZ+mub9h+PmRDEcKquaxcO0nMG84LBzhp0F3yTB/OPz4CfR62U+XHgZ92tdj9IKfM+3SFxtj9G5XLyzXkehQS5SIiIiIlB6xZeGM/4OBH0LVBsHyX5bDS51hxbiwTIXeqGY8w3q1JNYs7eXNGN6rJY1qxhf4GhI9aokSERERycKah8+JdggSKfVPgRsWwqyh8NUbvuzgHphxK/wwCy58DuJrFugSvdvVo1mdylz4/KJDZdOHdKRV3aoFOq9En1qiRERERLJQMa5MnjYpZsofBj1HQu/xUKFasHztB/BCB/hhdoEv0aBGfLb3pXhSEiUiIiIipVvzi/xU6Ed3DZbt+QPeuswv0puY/Ux7UvooiRIREREROewIuGoadH8UYssFy5e9AqNOg19WRi82KXKURImIiIiIAMTEwCk3wPXzoU7LYHnCOhhzFnz2BKQkRy8+KTKURImIiIiIhDr8OLh2DnS6DQjMrpeSBJ/+B147D7ZviGZ0UgQoiRIRERERSa9MOTjrYRgwE6qErOm0eSm8eCp89VZYpkKX4klJlIiIiIhIVhqe6qdCb9k7WJa4C965ESb3h71/RS82iRolUSIiIiKSJ3sTkzJs2e0v9ipUhUtGwyVjoHyVYPmad/1U6D/OiVpoEh1a0EBERERE8qT5Ax9lu7/DsE/T3N8w/PxIhlN4Wl7qF+mdfgNsWODLdv8Gb/SCk2+AMx+EshWiGqIUDrVEiYiIiIjkVpW60O89OPsRiI0Lln/+ErzcBX79OmqhSeFRS5SIiIiI5Mmah8+JdgjRFRMDHW+Bo7vA1MHw53e+/M/vYXQ3OOP/oMPNEBMb1TAlctQSJSIiIiJ5UjGuTJ62EqtOS7huHpwyJFiWchA+fgDGXQg7NkctNIksJVEiIiIiIvlVtjx0/x9c/Q5UPiJYvnEhvNiJsmumRC00iRwlUSIiIiIiBdW4K9y4GJr3DJYd2EnFGTfwTNnnOYzdUQtNwk9JlIiIiIhIOFSsDpeNhYtHQVzlQ8UXxS7mg3L308Z+jF5sElZKokREREREwsUMWveFGxdB/Y6HiuvaNibFPUTcilfAuSgGKOGgJEpEREREJNyqNYABM9l/2r9Icv4jd5wlU+GT+2DKIDiwK8oBSkEoiRIRERERiYSYWA50uJ3LE//F765qsPzbafByV/h9TdRCk4JREiUiIiIiEkHLXDPOPzCMxcnNg4UJ6/yaUqsmRi8wyTclUSIiIiIiEbaNKlx18H72d7gzWJi0D6ZfDzNug4P7oxec5JmSKBERERGRQpBCDAdOux+umAzlqwZ3rBgLY86Cv9ZHKzTJIyVRIiIiIiKFqcnZcMMCOPLEYNlvX8Oo0+H7D6IXl+SakigRERERkTDZm5iUYct0f8Uj2XvVTDjpuuDOAzth4uXw8QOQnIQUXWWiHYCIiIiISEnR/IGPst3fYdinae5vGP441DsZ3rsVDu7xhYuegc3L4NJX4bAjIhWqFIBaokREREREoqnlpXDdXKjVLFi2aTGM6gzrP4teXJIltUSJiIiIiITJmofPyd+BtZrC4E9h5h3w9du+bM+fMP4i6PpPOPVOiFH7R1GhJEpEREREJEwqxhXg43VcPFw8Cup3gA/vgeREcCnw6X9g8+d+X8Xq4QtW8k3prIiIiIhIUWEG7QbCNbOhaoNg+brZMOo02LIierHJIUqiRERERESKmiNPgOvnQ9PzgmU7N8Or58AXo8G56MUmSqJERERERIqkCtWg71tw5kNgsb4s5SB8cDdMvQYO7IpufKWYkigRERERkaLKDE69HfrPgEq1g+XfTIXR3eCP76IWWmmmJKqEMrPTzOxdM9toZs7MHox2TCIiIiKSTw07wfULoGHnYNm2H3witert6MVVSimJKrkqAWuAe4DfohyLiIiIiBRU5dpw9TvQ+a5g2cG9MP06mHE7HNwfrchKnVKTRJlZdTP7I9Aq830U42hrZveZ2TQz+yUQT65+482svJk9ZGY/mNl+M9tqZq+aWd30dZ1zHzjnhjrn3gYOhP2BiIiIiEjhiy0DZzwAV0yC8lWD5Steg1fPhu0bohVZqVJqkihgBFAz2kEA/wcMAy4GjsztQWZWHpgDPIBvZXoX2AwMBFaaWePwhyoiIiIiRVKTc+D6z/wsfql+XeWnQf/+g+jFVUqUiiTKzM4A+gOjox0LsAR4GOgB1MnDcfcDHQPHN3HO9XHOnQzcBdQCXg13oCIiIiJShFVrAIM+gvbXBsv274SJl8PHD0ByUvRiK+FKfBJlZhWAl/Djg57Ix/GdzOzSHOpUMbMHzCzH59M596hz7t/OuZnOud9zGUNZ4JbA3SHOud0h5xsBfA2cZmZtc3M+ERERESkhypSD85+ES8ZA2fhg+aJnYPyFsEtD4yOhxCdRwL+BxsCNwMG8HBhIXt4CJphZzyzqVAI+BB4CripQpFk7FagK/OSc+zKT/VMCtz0idH0RERERKcpaXgrXzYWaTYNlGxfBS51h/WfRi6uEKtFJlJm1wnd3e805l+ffHufcQaAvsB9428zOC91vZhWBD4AOwHjg9QIHnbnWgduVWexfma6eiIiIiJQ2tZrC4E+hZe9g2Z4/YPxFsOBJSEmJXmwlTIlNogJd60YDO/DTfOeLc24JcAGQBEwzs7MC5y8PzAA6AxOBQc45V8Cws1I/cLsli/1b0tXDzCqZWRszawPEAXUC95vndDEzq2pmDc2sIVA2RX9wIiIiIsVDuUrQ62U4fwTExvkylwJzHoYJfWHvX9GNr4QosUkUfgzRScA/nHMJBTmRc24+0DNw910z6w68A3QDpgNXO+eSC3KNHFQK3O7NYv+edPUA2gFfBrYjgOsDP+dmupbbgfWB7diEhAI9fSIiIiJSmMyg/TVwzWyoWj9Yvu4jGHU6/LIierGVECUyiTKzesAjwHzn3NhwnNM59zFwCVAGPwbqHGAm0Mc5F+mpTyw1jBz2H+Kcm+ecs0y2hrm43tNAo8C2rkaNGvmJWURERESi6cgT/DToTboHy3Zugle7wxejIWKdqEq+EplEAS/gu7DdGObzfgR8HnL/0cC4qUjbFbiNz2J/xcDt7iz254lzbodzboNzbgNwMCampP6aiIiIiJRwFapB3wlw5oOQOpF0ciJ8cDdMvRYOhOXjY6lTUj8dX4Dv+vaimc1L3fBjlwDqh5RXyvIsIcwsFj9xxKnAKiAZeC8w5ijSNgVu62axv266eiIiIiIiXkwMnHoH9J8B8YcHy7+ZAqO7wR/fRy+2YqqkJlHgpwQ/Pd12cmBfhZCyMjmdKDBJxav4mfoWAp2Aq4EqwMdmdnyYY09vVeD2xCz2p5Z/HeE4RERERKS4angq3LAAGpwaLNu2FkZ3ha8nRy+uYqhEJlFZjAUy/BgfgLUh5TuyO5eZGX6x3n7AUuA859we59wEYBBQA5hjZsdF7hGxCNgJNDazEzLZn7oY8MwIxiAiIiIixV3lOtDvXTj1zmDZwb0w7VqYeSckHYhebMVIiUyiwuw5YDCwAujunEsdn4Rzbhx+1rta+ETq2EgE4JxLBJ4P3H3ezA6NjTKzO4FWwELn3LJIXF9ERERESpDYMnDmv+Hyt6F8lWD58jEw5mzYviFqoRUXSqKyYWZlgMb47nRnO+d2pq/jnBuNn069ClAnF+c838yWpm6B4rjQMjM7P5NDH8FPatERWGdmbweOfxJIAAbm5zGKiIiISCnVtLufve+INsGyX7+CUafB2g+jFVWxoCQqG4Gpy3sBZzjnslyZzDk3EmjmnFuQi9PWwo/NSt3AT1EeWlYrk2vsB7oC/8FPmtETaAiMA05wzv2YqwclIiIiIpKqWkO/nlT7a4Nl+3f6hXk/eRCSI72ST/FkTvPDSzbM7NvmzZs3//bbb6MdioiIiIhE0teTYcatfoxUqganwqVj/FiqEqZFixasWbNmjXOuRV6PVUuUiIiIiIhAq8tg8Fyo2TRYtnEhvNQZ1uemw1XpoSRKRERERES8w5vB4E/h+EuDZXv+gPEXwoIRkJKSr9PuTUzK01bU5bhGkoiIiIiIlCLlKsElr0CDDjBrKCQngkuBOQ/B5s+h54tQsXqeTtn8gY/yVH/D8MzmWSs61BIlIiIiIiJpmfnJJgZ9BFXqB8t/mAUvnw6/rIxebEWAWqJERERERCRzR50I18+H6TfAukBr0o5N8Oo50H0YtLvGJ1w5WPPwOWnu79x3kA7DPj10f8nQblSpUDasoUeSWqJERERERCRrFavD5RPhjH+DBdKH5ER4/y6YNhgO7M75FHFlMmzZ7S/qlESJiIiIiEj2YmKg853Q7z2IPzxYvnoyjO4Gf66NXmxRoCRKRERERERyp1FnuGEBNOgULNu2Fl7u6teZKiWURImIiIiISO5VruNbpDrdHiw7uAemXeu7+CUdiFpohUVJlIiIiIiI5E1sGTjrIT9WqnyVYPmyV/ykE9s3Ri+2QqAkSkRERERE8qfpuXD9Z3BE62DZ1i9h1Gmwdlb04oowJVEiIiIiIpJ/1RrCoNnQblCwbP8OmNAHPnkIkpOiFVnEKIkSEREREZGCKVseLngKeo2GshWD5QtHwOs9YdfvUQstEpREiYiIiIhIeLTqDYPnQs0mwbINC2BUZ9iwMHpxhZmSKBERERERCZ/Dm/lE6vhLgmW7f4dxPWDhU5CSEr3YwkRJlIiIiIiIhFe5SnDJGDjvCYgp68tcCnzyIEy8wo+ZKsaURImIiIiISPiZwUmD4ZqPoEq9YPkPH1J5bDeOt5+jF1sBKYkSEREREZHIOaqtnwb92LMPFcXs3MTUuAe5InYO4KIXWz4piRIRERERkciqWB0ufxvOeADMpyDlLIn/lR3DU2VfgMQ9UQ4wb5REiYiIiIhI5MXEQOe7oN+7pFSsdai4pa2nuLVGKYkSEREREZHC0+g0dv9/e/ceLVdZJmj8eQkESKJcAt3AAEJH2gEa5KY0dy8sYBBaRBpoVBoYL4um6UG8TIPDTewWnOWlFdBpXAiKIA4DchGWTWMHQUERHKQblBgIBhGZBBKuSSB5549vF1SKc07VPufUqVN1nt9ae+2z9/72t9/aO2ul3vou+7gf8tNV/5kXcm1OfOkUmD6r11HVsmavA5AkSZI0teTrNuWYFZ9i23iUebl5r8OpzZYoSZIkSRNuJdP49/yTXocxKiZRkiRJklSDSZQkSZIk1WASJUmSJEk1mERJkiRJUg0mUZIkSZJUg0mUJEmSpAn16OLnR9ye7EyiJEmSJE2Y7/58Ie+58Cer7XvPhT/huz9f2KOI6jOJkiRJkjQhHln0PKddcz8rM1fbvzKT0665n0cW9UeLlEmUJEmSpAlx1d0LWbkqhzy2clX2TWuUSZQkSZKkCfHY0y+0Of7iBEUyNiZRkiRJkibE5hvMaHN83QmKZGxMoiRJkiRNiKPesgXT1oghj01bIzhyty0mOKLRMYmSJEmSNCG23mgmnz18B6bF6onUtAjOO3wHtt5oZo8iq8ckSpIkSdKEOXK3Lbj2pD1X23ftSXvyl33SCgUmUZIkSZIm2Btmzxxxe7IziZIkSZKkGkyiJEmSJKkGkyhJkiRJqsEkSpIkSZJqMImSJEmSpBpMoiRJkiSpBpOoARUR+0bEdRHxaERkRJzd65gkSZKkQWASNbhmAQ8AnwSe6HEskiRJ0sAY2CQqIk6NiGsiYl5ELI2I5VWrzGURsX0P49o1Iv6+iu13VSvRsg7PXScizomIhyJiWUQ8HhGXRMTmrWUz86bMPC0zrwKWj/sHkSRJkqaoNXsdQBedDswEfgncX+3bHjgWODoiDsvMm3sQ1xnAu+ueFBHrALcCewK/B64DtgKOBw6JiD0yc/44xilJkiRpCIOcRL0buCczV2vliYgTgYuAr0fElpm5coLjuhO4D7i7Wjrtanc6JYG6EzggM5+D0uIGfB64BNhv3KOVJEmStJqB7c6XmT9uTaCq/V8FfgNsBrypXT0RsVdEHNGmzHoRcWZEtL2fmXl+Zp6VmTdm5h/ala/qXws4udo8qZFAVfV9gdLatm9E7NpJfZIkSZJGb2CTqDYarU8rRipUJS9XAFdGxGHDlJkF3AycA7x/HGNstjewPjA/M38xxPGrq/WhXbq+JEmSpMqUS6Ii4lhKC9RDwMMjlc3Ml4CjgWXAVRFxcEtdM4CbgD2AbwLf6kbMwJur9b3DHL+3pZwkSZKkLhnkMVEARMQnKBNKzAS2rf5+HDgmM1e1Oz8z74yIQyjJ0jURcWhm3lJN9HADsA/wHeCEzMwufYwtq/Vjwxx/rKVco4XsjdXmdGCTiNgJWJGZD4x0sYhYn9LyBbDWqlVtb5MkSZI0ZQx8EgUcCLyzaXsh8IHMvKfTCjLztqo73w3AdRFxOHAK8A7g2qq+bk5QMatavzDM8edbygHsBvxb0/ZHquVRyqx+IzkFOKuxsXjx4g7DlCRJkgbfwHfny8z9MzOADYB9gV8DcyPiUzXruQV4LyXxvJmSnN0IHJWZL49v1K8RjTDaHH9FZs7NzBhi2aqD630J2Lpa5s2ePXs0MUuSJEkDaeCTqIbMXJKZtwMHA/cA50bEW2pW8wPgp03b51fjprrt2Wo9c5jjM6r1c8Mcr6W6VwsycwHw0hprTJl/JpIkSVJbU+7bcZX0XEVpvel4NruImEaZOGJvynueVgLXV+OMuu231XrzYY5v3lJOkiRJUpdMuSSqsqhab9xJ4er9T5dQZuq7A9gL+ACwHnBLRPxZN4Jscl+13mWY4439v+xyHJIkSdKUN1WTqP2q9fx2BSMigK8BxwJ3AQdn5vOZeSVwAjAbuDUitu1WsMCPgaXAnIjYeYjjjZcB39jFGCRJkiQxoElUROwTEUdFxJot+9eKiJMprUgvUrr1tfMV4EOUcVQHZWZjfBKZeRllxruNKYnUNuP1GZpl5grggmrzgoh4ZWxURJwK7AjckZl3d+P6kiRJkl41qFOczwG+ASyKiHuAxcBGwA7AppSX5x6XmQtHqqRKwuZQutMdkJlLW8tk5sURMR34HLAJMK9Nne8CzmjZPT0i7mraPjczv99S5jPA/sCewLyIuB14A7B79fmOH+m6kiRJksbHoCZRtwH/SOm2tyMlgVoBLACuBr6cmb9pV0lmvly9E2pGZj41QrkLI+L6dklZZWNK4tMsWva9ZqxWZi6LiLcDpwHHAIcBTwOXAWd0eG1JkiRJYzSQSVRmPgLUeg/UCHW9SOn6165cR0lMZl4KXDqGWM6sFkmSJEk9MJBjoiRJkiSpWwayJUqSJEnS5PHCipfbbq81LV7ZnjF9cqcpkzs6SZIkSX1vuzN/MOLxPT77w9W2F5z3rm6GM2Z255MkSZKkGmyJkiRJktRVD3z6wF6HMK5MoiRJkiR11WQf41SX3fkkSZIkqQaTKEmSJEmqwSRKkiRJkmowiZIkSZKkGkyiJEmSJKkGkyhJkiRJqsEkSpIkSZJqMImSJEmSpBpMoiRJkiSpBpMoSZIkSarBJEqSJEmSajCJkiRJkqQaTKIkSZIkqQaTKEmSJEmqwSRKkiRJkmowiZIkSZKkGkyiJEmSJKmGyMxex6BJLCKeWXvttV83Z86cXociSZIkjZv58+ezfPnyZzPz9XXPNYnSiCLiCWAGsLDXsWjcrQHMBhYDq3ocyyDy/q5u0O9HP3++foh9ssXY63gav2zO78G1pUGyBfBCZm5S90STKGmKioitgEeArTNzQW+jGTze39UN+v3o58/XD7FPthh7HU9E/AdAZm4/0deWVDgmSpIkSZJqMImSJEmSpBpMoqSpawlwTrXW+FuC97fZEgb7fiyhfz/fEiZ/7EuYXDEuYXLFI2mCOSZKkiSpjzgmSuo9W6IkSZIkqQZboiRJkiSpBluiJEmSJKkGkyhJkiRJqsEkSpIkSZJqMImSJEmSpBpMoiRNChGxb0RcFxGPRkRGxNm9jmmQeH+nFp+3JHWXSZSkyWIW8ADwSeCJHscyiLy/U4vPW5K6yCRK6pKIODUiromIeRGxNCKWV78KXxYRPXtBYkTsGhF/X8X2u+pX6mUdnrtORJwTEQ9FxLKIeDwiLomIzccaV2belJmnZeZVwPK650fEhhHxZPV5fjXWeEZrUO/vZBERm0TEF6t79GJEPBUR90TE53oUj89bkqagNXsdgDTATgdmAr8E7q/2bQ8cCxwdEYdl5s09iOsM4N11T4qIdYBbgT2B3wPXAVsBxwOHRMQemTl/HOOs6wvARj28fsOg3t+ei4g9gJuA9SmtLNcDrwO2A06ltLpMNJ+3JE1BJlFS97wbuCczV/tVOiJOBC4Cvh4RW2bmygmO607gPuDuaum0q8/plC98dwIHZOZzUFrcgM8DlwD7NQpHxIbAhm3qfCYzn6wV/RAi4p3AXwP/DHx4rPWN0cDd38kgIjajJFBrA4dn5rUtx9/ak8B83uoDEXEA8FnKD3mLgW8AZ/Xg/x9pcGSmi4vLBC/APCCB7ToouxdwRJsy6wFnAmuMIpYElrUpsxbwdFV25yGO31cd27Vp39nVvpGWS4e53gLg7A7jX7e6n/8BbFPV+6san9/72wcL8M3qM/3tGOvxebtMqQXYGVgBfBHYFngvsAQ4r9exubj08+KYKKk3Gr/+rRipUESsBVwBXBkRhw1TZhZwM3AO8P5xjLHZ3pQuVPMz8xdDHL+6Wh/atO/TlC+LIy0njENsZwFzgBOBl+qc6P3tDxGxAXAksBT4+hjq8XlrUhntmLqa4+k+DtyfmR/NzAcz8/9QfiT4u+rfu6RRMImSJlhEHAu8CXgIeHikspn5EnA0sAy4KiIObqlrBqWL0x6UX+q/1Y2YgTdX63uHOX5vSzkyc1VmvtxmWTWWoCJiR+BjwDcy80d1z/f+9o29KN347gBeiogjIuJLEXFhRJwcEX/cSSU+b01CZ1C62b0H2KyTE5rG051JmYXxOmAhZTzdvRExp+WUvSj/rpt9n9KKv+uoI5emOMdESV0WEZ+g9EOfSelKsT3wOHBMJ196MvPOiDiE8p/gNRFxaGbeUv1HegOwD/Ad4ITMzC59jC2r9WPDHH+spVxt1S+ib6w2pwObRMROwIrMfGCI8msAF1O6pYx6QgHv79D3d5JpzGb5B+B2SpLT7LMRcXxm/u92Ffm8++J5TyWjGVNXazwdsCll8pJmjet0lLhJei1boqTuO5Ay6cERlC+DCykJ1D2dVpCZtwGHVZvXRcRBwPeAdwDXAh/I7g4QbnT5eGGY48+3lBuN3YBfVMumwEeqv1t/QW04GXgr8InMXDyG63p/J78NqvWxwI7AfwU2BramzMo4E7i8aplsy+etySIzz8/MszLzxsz8Q7vyVZfUk6vNkxoJVFXXFyizwe4bEe1amLJlLakmkyipyzJz/8wMyhfBfYFfA3Mj4lM167mFMiB4TcqYjQOBG4GjMvPl8Y36NaIRRpvjo5aZczMzhli2es3FIrYAPgPclpmXjvXa1fW9v5PXtGq9JnBqZl6SmYsyc0Fmfowyhmg6NVokfd7qU6MZT/d7SiLdbNOmY5JGwSRKmiCZuSQzbwcOBu4Bzo2It9Ss5gfAT5u2z6/GeXTbs9V65jDHZ1Tr54Y5Pt4uonxpPnGc6/X+Tk6N+7MKuGyI45dU67fVrNfnrX5Tezwd8GPgv7SUOxh4kfJ/kaRRMImSJlj1Je0qyq/Nh7Yp/oqImEYZ6L43pQ/9SuD6apxDt/22Wg8181Pz/t8Oc3y8HULpCvXViJjbWChjWQC2bNrfUZco7++ktqBaP5GZy0c4/kedVujzVp8azXi6zwM7RsTnI2LbiDicMtvjl5u7A0qqxyRK6o1F1XrjTgpXkyhcQplZ7A7KbEsfoLzP5paI+LNuBNnkvmq9yzDHG/t/2eU4mq1PGTzdvOxeHVu3aV/bCXS8v5Neo9vSBhExVFe32dW6oy+EPm/1sdrj6TLzXsoPdm8H/i9wAaU1v1aXckmrM4mSeqMxc9L8dgWrL41fowyqvws4ODOfz8wrKe+FmQ3cGhHbditYSneQpcCciNh5iONHVOsbuxjDK4YZ6xGUiQYAft20f8lIdXl/J7/MvB94hJIc7z5EkbdV6+G6OL3C560+N6rxdJn5g8zcJTPXzszNMvNTXZ48RRp4JlFSF0TEPhFxVESs2bJ/rYg4mfKr94uUbn3tfAX4EKXv+kGZ2RhPQWZeRpl1a2PKF79txuszNMvMFZRfLwEuiIhXxnJU0+ruCNyRmXd34/pd5v3tD+dX6y9HxEaNndUsZB+rNr/WQT0+b/Uzx9NJk4TviZK6Yw7wDWBRRNwDLAY2AnagzIq0DDguMxeOVEmVhM2hdP85IDOXtpbJzIsjYjrwOWATYF6bOt9FecFjs+kRcVfT9rmZ+f2WMp8B9qe8n2ReRNwOvIHSMrCY8qLHvuL97SsXA+8E/hL4dUT8hNJlaU/KJCMXZ+bVI5zv89YgcDydNEmYREndcRvwj5RueztSEqgVlAHwV1MG9P6mXSWZ+XI1CHhGZj41QrkLI+L6dklZZWNe2yUqWva9ZqxWZi6LiLcDpwHHUN6z8zRltrQzOrz2pOL97R+ZuSoijgbmAh+kvNMpgZ8DX8vMb3VQh89b/c7xdNIkEd17IbskSZI6FREJLM/MdYY5Ph14kjIJyi6t74qKiPsoP9y91e6gUnc5JkqSJKkPOJ5OmjxsiZIkSeqBIcbU7U7ppvqzpn2rjamLiHUo3Vp3B34PtI6n+/NOuotLGhvHREmSJPVG7TF1jqeTJgdboiRJkiSpBsdESZIkSVINJlGSJEmSVINJlCRJkiTVYBIlSZIkSTWYREmSJElSDSZRkiRJklSDSZQkSZIk1WASJUmSJEk1mERJkiRJUg0mUZIkSZJUg0mUJEmSJNVgEiVJkiRJNZhESZIkSVINJlGSJEmSVINJlCRJkiTVYBIlSVKPRMS+EXFdRDwaERkRZ/c6JklSeyZRkiT1zizgAeCTwBM9jkWS1CGTKElS34uIC6qWnObl2Yi4LyJOj4h1W8q/bYjy6zcd36raN7ebcWfmTZl5WmZeBSwfqWxLrHuMUO7IpnILxjvmdiJiSUusx010DJLUbWv2OgBJksbBm6v1TcD/q/7eAtgP+AfgXRGxT2auajlvPnBH9feKrkc5ft4H3DnMsfdPZCBDuAKYAezEq89FkgaKSZQkqa9FRAA7VpsfzszfNR3bH/gXYE/gQODmltPvyMzjxjGWDYEN2xR7JjOfHOUlllMSv6Mi4pTMfLnl+rOBg4B7gV1GeY0xycy/qWI5G5MoSQPK7nySpH63NfB6SnLyu+YDmfmvwM+qzR1bT+yCvwPmtVk+N8ZrfBvYiJIUtjoKWAu4fIzXkCSNwCRKktTvGq0dDw5zvNG9b8loKo+IdSPivGoGveUR8ZuI+O9VC1irT1OSmJGWE0YTR5NvA8nQ3fbeDzwHXDfMZ3llrFdEvD4i/ikiFkbEsoh4MCI+GhFDfjeIiC2rsWfzqvKLI+JnQ405k6RBZ3c+SVK/a5dEvala3z2KuqdTugNuT2nRepAyzuo84HXA/2guXI25ah13Na4y89GI+DHwFxExKzOfA4iIrYE9gG8CL7SpZm3gh8Ccaj0deCfwBUqL3fHNhSNiX+B6YD3gYUqSNhPYjjLm7ApgwTh8PEnqC7ZESZL6XSOJeqD1QES8F9gGuD4z7x1F3XtQWn3+NDMPysyDgH2Al4GPRsSsUcbciG9WROwUETtREplNqu3t2px6OWXyhsOb9jVapr7dwaX/nNIqtk1mvjczDwV2AB4HjouIv2iKcQPgakoC9VHgjZl5VGYekpl/Qkkqn+7gmpI0MEyiJEn9bqdq/UpLVERsGBEfBC4BvgccM8q6VwEfzMxFjR2Z+XPKBBUzgN1GWW/DbsAvqmVT4CPV3ze1Oe+7lNkE39e0732Ud03d2uG1P97yueYD51abJzWV+xCwMXBjZn4pM7O5ksz8UWYu7fCakjQQ7M4nSepbEfF6YKtq84Yhhin9U2aeMoZLLMjMh4bY39i36RjqJjPnAkONrWp33tMRcRNwaERsQpnO/U3AFzNz5dDDtVbzVGbeMsT+K4CvAntGRFQJ0/7Vsf9VN05JGlQmUZKkftboyreU0uIEZazOW4Etgf8WEXMz83uvPbUjjw2z/7lqvfYo6x0PlwOHAUdTZihs7OvEo0PtzMxnImIJsD5lxsOllAQNytTqkiRMoiRJ/a2RRN3e/L6niFiL0iVuf+BkXk2w6sr2RXrmRsqMg8cCmwEPjnLcV6vhmrEm872QpAnlmChJUj/bqVqvljxk5kvAadXmfhGx0UQGNREyczllwoedgT+m3ruhthxqZ9U9cj3geeCZavfCav3G0UUqSYPHJEqS1M8aLVGvaYGpJoD4FTANOHgig5pA3wQWA4vobFa+htkRsf8Q+/+qWv+kaQKJf63WHx5diJI0eEyiJEl9KSKmUd7fBEMkUZXGS2cP6X5EEy8zb8/MjTJz48wccpzTCP5nRMxubFTvmTqj2ryoqdzXKUnaoRHxt60vGY6IfSJivdHEL0n9yiRKktSv/hRYF1iUmQuHKXN9tT6gGiel4i7K9O3zIuLqiLge+HfgPwGXN0/EkZlPAUcCzwJfqc65KiJuiIiHgR8BG0z0B5CkXjKJkiT1q2G78jW5C3iSMs5nn65H1D+WA+8ArqS8UPhAytinjwPHtRbOzH+jjD/7Z8qkVIdRXtj7JGXs2RPdD1mSJg9n55Mk9aXM/A7wnTZlVlEmXahb9wJGeH9TZp4NnF233rHIzI7fJ5WZT9Dm/VPVC3JPYvUX645U/mHKy4AlacoziZIkTWV7R8Sl1d9/k5kv9DKYQRARFwEzeHXmREkaOCZRkqSpbE61AJwCmESN3TGU7pOSNLBMoiRJU05mzqVNdzeNTmau3+sYJKnb4tXXQEiSJEmS2nF2PkmSJEmqwSRKkiRJkmowiZIkSZKkGkyiJEmSJKkGkyhJkiRJqsEkSpIkSZJqMImSJEmSpBpMoiRJkiSpBpMoSZIkSarBJEqSJEmSajCJkiRJkqQaTKIkSZIkqQaTKEmSJEmqwSRKkiRJkmowiZIkSZKkGkyiJEmSJKmG/w+FestvmKimvgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure(dpi=150)\n", "# curve_fit requires initial guess (p0) to get resonable fits, sigma takes in the y errobars \n", "popt, pcov = curve_fit(model, x, y, p0=[12,0.5], sigma=yerr)\n", "\n", "log_M200m, c = popt\n", "sig_logMh, sig_c = np.sqrt(np.diag(pcov)) # pcov is the covariance between the parameters\n", "\n", "# rough estimate of chisq\n", "chisq = np.sum((y-model(x, log_M200m, c))**2*1.0/yerr**2)\n", "dof = 10 - 2 # 10 datapoints - 2 parameters\n", "plt.errorbar(x, y, yerr=yerr, fmt='.', capsize=3, label='Data')\n", "plt.errorbar(x, model(x, log_M200m, c) , fmt='-', label=r'Fit, $\\chi^2 /{\\rm dof} = %2.2f/%d$'%(chisq,dof))\n", "\n", "plt.title(r'$M_{\\rm 200m} = %2.2f \\pm %2.2f \\times 10^{14} \\,h^{-1} M_\\odot,\\, c = %2.2f \\pm %2.2f$'%(log_M200m, sig_logMh, c, sig_c))\n", "plt.legend()\n", "plt.xlabel(r'$R[{\\rm h^{-1}Mpc}]$')\n", "plt.ylabel(r'$\\Delta\\Sigma [{\\rm h M_\\odot pc^{-2}}]$')\n", "plt.xscale('log')\n", "plt.yscale('log')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When we compare our halo masses with the ones using SDSS shape catalog [arXiv:1707.01907](https://arxiv.org/abs/1707.01907), we found that within the error bars our results are in agreement with the reported ones as given in subpanels of fig 4 and fig 5 in [arXiv:1707.01907](https://arxiv.org/abs/1707.01907).\n", "\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3.8 (py38)", "language": "python", "name": "py38" }, "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.8.5" } }, "nbformat": 4, "nbformat_minor": 4 }