Skip to content

Instantly share code, notes, and snippets.

@smsharma
Last active March 25, 2019 16:23
Show Gist options
  • Save smsharma/d399da0d8ad905f9f431553e34c33e3c to your computer and use it in GitHub Desktop.
Save smsharma/d399da0d8ad905f9f431553e34c33e3c to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this notebook we will make heatmaps to see which galaxies (by redshift and virial mass) are most important for DM studies."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import os, sys\n",
"from itertools import compress\n",
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"import healpy as hp\n",
"import matplotlib.pyplot as plt\n",
"import astropy.units as u\n",
"from astropy.cosmology import Planck13, z_at_value\n",
"from astropy.coordinates import SkyCoord, Distance\n",
"import scipy.interpolate as ip\n",
"from scipy.interpolate import UnivariateSpline, interp1d\n",
"from tqdm import *\n",
"import matplotlib.gridspec as grd\n",
"\n",
"from units import *\n",
"import cart2sph as cs\n",
"\n",
"%matplotlib inline\n",
"%load_ext autoreload\n",
"%autoreload 2"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Define quantities we want to plot\n",
"mvirvals = #\n",
"zvals = #\n",
"Jfunc = # Function to get Jfactor from M and z"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# z and M bins to use\n",
"z_bins = np.linspace(0,.06,100)\n",
"M_bins = np.logspace(12,16,100)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"J_heatmap = np.zeros((len(z_bins),len(M_bins)))\n",
"\n",
"for i in tqdm(range(len(mvirvals))): # Loop over all halos\n",
" Jf = Jfunc(mvirvals[i],zvals[i])\n",
" for j in range(len(M_bins)-1):\n",
" if (mvirvals[i]/M_s >= M_bins[j]) and (mvirvals[i]/M_s < M_bins[j+1]):\n",
" for k in range(len(z_bins)-1):\n",
" if (zvals[i] >= z_bins[k]) and (zvals[i] < z_bins[k+1]): # If passes selection, add to bin\n",
" J_heatmap[j,k] += Jf"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"x = z_bins\n",
"y = np.log10(M_bins)\n",
"X, Y = np.meshgrid(x, y)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAE7CAYAAADQCfvFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXl8XHd57/9+Zkab5UWW92y2R86eQCIrSxMgJJYSNrMk\ndlKWQltiKfRSyi3FwrS9r/6gxZHpetvLL5JpCy1w61iBBAElkQylEAjYUjYSskneEjuxLVmytWtm\nnvvHOTMejUbSjGbRjPS8X6/zkuac7/me75ljz3z0rKKqGIZhGIZhzCU8s70AwzAMwzCMdGMCxzAM\nwzCMOYcJHMMwDMMw5hwmcAzDMAzDmHOYwDEMwzAMY85hAscwDMMwjDmHb7YXkC2WL1+u69atm+1l\nZIz29nY2btyY0Ni+vj6WLFmSN2NPnjzJihUrEhqbyXXkwv0lM2+y4zM19siRI1x00UWzuoZMzp2p\n55cr70W+3V+y70V7e3u/qi5K+IRJuE5E+1KdZAa8DI+q6jtm4dK5j6rOi23jxo06l/F4PAmP3bZt\nW16NTfbZ5cKaM3V/ycyb7PhMjV2+fPmsryGTc2fq+eXKe5Fv95fsewEMaBq+Yy4Gbc3i9gXQd4EC\nL6dj/XNxmzcWHOMcmzdvzquxyZILa87U/SU7by7cXzJ/TWfyPc6355cr70Wm5s2FselEyK5L5K3u\n9gOYDcNRXiCOgJ37VFVV6YEDB2Z7GRnD6/USDAZnexkZoaqqirn87Oz+8hu7v/xGRAZVtTTVeS4V\n0cZ0LChJboV2Va2ahUvnPGbBmSP4fHP3UdbW1s72EjKK3V9+Y/eX95xMxyQCFKRjogwhItVAjarW\nxzkW/ZDLVHXXNHP5gTqg093VpaptaVtsmpg3FpyLL75Yb731VjZv3jxrJsxM0NLSQktLC1/72tf4\n3d/93Tl3f4ZhGJkg/Nm5e/fuV1T14lTnu0xE/zkdC0uSt0xjwQkLG6ASR4jUxRxviBY9IlIJVE8m\nclxx06iqNe7r7TjCqSb1u0kv80bgzHUXVWlpKQMDA7O9DMMwjLxCRNLi4rlcRL+WhvUky40JuqhE\npAHHOhMrcFpjxYmI7FXVrZPMsxfYo6rN7usyAFXtnek9ZAqrg2MYhmEY8xgR2RsWKq7FZ88k48qA\nLUDEHaWqvbkobsBicAzDMAwjZXI9BmcK6oB24KCI7MRxYzVPMjZsKfK7riqA6+LF9eQCZsExDMMw\njHmKqnYB9ThWmQbgnimGl7k/y1W12RVCnSIyGwlk02ICxzAMwzDmKa44aXNjbmqAajfOJh69ADEZ\nUweAnEy1MxeVYRiGYaRItgr9NQMPjd+1fKZzuRlTva4VB1VtE5H1OO6qsjixNT1xpul154o3flYx\ngWMYhmEYKZKtGJwPuluYK+BUCtOVc66WDeAEDYvIg/EGq2qHiPTGiBnLojIMwzAMI6c4AIxLB3cz\npSKZUSJSGVMIsAm4O+p1NTBlYcDZwiw4hmEYhjFHCRfuw0nvDhfma1PVDtdaU+fWyOnGcTeVxWRF\nVeMEHjcBqGq9iDS48wAsy9UsqnkjcPr6+qitrZ1zlX7D1TjHxsbm5P0ZhmFkgvBnJ5B4N9gpyHaz\nzURR1Q6gg0msLFFZVJOdvyv23FwVNLFYJeM5glUyNgzDSJ50VTK+WkS/nY4FJckl1mxzUnJRcBpG\nyrS1teH3+/H7/ZF9GzdupKpq/OdARUUF27dvH7evrq6Onp4e/H4/vb29lJWV0dDQkJV1Z5Lm5mZ6\nenro7Oxk2bJlE+47U9TU1ES6UdfW1s6J99IwYsnjQn9zlpwTONN0PM2LDqbG7FFfX09vby9NTU3s\n3bt3nMDp6Oigo6Nj3PjW1tYJc/T09NDc3ExZWdmc+ULu6urC7/ezZcsWwBF71dXVVFZWZvzadXV1\n7N27l7KysukHG0aeYgIn98iZLCoRqXYDneo5Vy0x+ni4g2m9qja5Y/LCD2hkj4aGhriCpKuri4aG\nBlQVVaW9vZ2Ghgaqq6snjC0vL0dVOX369JwQN+CIu/r6c/9dqqqqaGvL3t8GJm4Mw8g2OWPBcS0x\nbeGOp3GGNADR5aCb3M0wpsXv91Nb62Q69vb20tjYSGPj5NXFm5qa6O3tpbu7mx07duT9F/SWLVvG\nibmuri5qamqmOCN9dHV10dzstLbp6emJPAfDMGZOq7uRpiDpuUjOWHCmIt86mOYTNTU11NTUUFFR\ngYggItTU1NDV1ZX1tXR1dUWERTx6e3upr69n69atbN26lbq6uqTWGRYp9fX11NXVTTqupqaG2tra\nSIzKxo0bk7gLh7a2NrZu3UpFRQUVFRVs3LiRmpoadu1ykhF27dpFU1Ny+rypqYmamprIc1q6dGnc\nZ9XW1sbGjRsRESoqKiL3Gr7/8PiwuyrTVFdXs2XLFrZs2UJra+us/NsyjGzgy+L2TuBvncv2Zf7O\n8pOcseBMQ151MM0nwjEou3btor6+Hr/fHzcuJVM0NzfT2NhIWVkZbW1t9Pb2UlVVNSE2pLe3l/Xr\n13P33Xezd6/TJqW+vp6NGzfS3t4+LtZmKsLxOVNZb6K/+CsqKujq6qKjoyPheJWtW7fS3NxMbW0t\n7e3t46w/bW1tkTmTdX/V1tZSW1tLXV0dTU1NlJeXx31W1dXV7Nu3j02bNtHe3j7heENDQ+Q9TIVd\nu3bR3d0d91hNTU3EYhT9vvn9fpqbm7MW4GwYxjwmHJOQKxuuKypm3xZAgeqofbWx46baNm7cqHOJ\nZw6N6DOHRiKvFyxYkNJ8DQ0NCqjf7592bGtrq1ZXV2tra2tK14zF7/croO3t7ROObdmyRQHt7OyM\n7Dt9+rQCWl1dPW5seP/evXsnzNPY2KhlZWWTrqG9vX3cdRobGyddUzyqq6sV0O3bt097jYaGhoTm\njKWzs1Pd/w9x71FVdfv27XGPNTQ06OnTp1VVIz8zyd69e7W2tnbcuhobGzN+XcNIFOCApuG7682g\n3b7sb+la/1zc8sJFRRo6mJ48eZKqqqrIlqx7INf41+8f5K+++njkfkZGRrJ2f9XV1RFXQ01NTVbe\ny3AMR7SlJmwZSSZYtr29nfLy8gn7w9lVZWVl+P3+yJjOzk7KysoSst40NzdH1jKVdaaysjJucHOi\n+P3+yPk7d+6cdC2xLqiwVSn8vmWjLpTf7x/nDuzo6Ejp3g0jVZqamsZ9VpJCs0ojt8kXF1XKHUxX\nrFiRlQ/0bLHmvDWsOW8N//FF555KS0uzfn9hl0lzc3MklicTroepYjbKysro7e2NuJB27dpFZ6dT\nRaCxsZH9+/ePExthARPLtm3b2L17N5WVldTX17Nz5056e3vp6uqK6+aJR9jtlcgX+GQxQGH3WVVV\nFT09PfT29tLQ0DBBYDU0NLBx48ZI6nv08aampgnzd3R0sHXruJYz095Xc3Mze/bsidQDCq8nGRdm\nZWUlzc3NdHV10dXVFXGDGsZsEf7cCiMiqTSrjJoHfLPxjRqYhWvmC7NtQordiOOicvefxumREX5d\n6Sw/sXnnmotq13f6dNd3+iKvs+mimoyw62r79u0zcn9M5qIKu3Tc5z2OsrIyBdLuLpsJ4fVv2bJl\nRufX1tZOcMOFn0v0vjCVlZVxXXR+vz9l91PYDRZNe3u7VlZWpjSvYeQapMnFc42gfcXZ39K1/rm4\n5YuLCvKog+l8Jey6uueee9i2bduUmUrGeMIZZJWVleMsHOG/NOMFRe/YsQNwXHRhK1fYNZVqWnt4\nvrq6uoj7rrKycs7UBTKMdCMCBd7sb8bk5IyLaqqOp5BfHUznO11dXfT29qbNFZHIl3W8uJps4/f7\nI66YWMLZYnBOPJSVlVFVVUVjY2MkdiecCh9NZWUly5YtmzDnli1bItdsaGigsbGRnTt3sm/fvpTv\nJVzluKmpKRJjVVlZye7du1Oe2zAMIxvkjMDRaTqeumNM0KSJ5uZm/H5/Wkv1h9sjbN26Na2p5lMJ\npXDNnFyI66irq6OtrW1COwggUgcGQEQAJ44mNmbJ7/cnZSUJ1/QJ18jx+/1pK0rY3t5OW1sbra2t\ndHR00NbWxqZNmzh48GDeFz40DGPuk08uKiNNNDc3s3Xr1rgF7GbyxbVr1y5qamoidVkyUak2LMTi\nFQFM55d6KkRXC07WPRc+L16geG9vbySLLJborKitW7dG3Fap0tbWxq5du6iurqahoYHW1lba29sj\ngdeGYYwnHGScre2HCn84Clgl40kxgTMPCX8hRltvwhaXe+65J6E5oqsKV1ZW0tramtHKuGHXSHRK\nevj3qYr2ZZvW1laqq6tpampi69atcasMx8Pv97N9+/ZIIcJotm3bNqWACwvKysrKtFnkent7I5lk\n0WuM/mkYxjkEKPBlb/tACXx1MTAHKhmLyGIRuVNE7heRPSLyqLvtEZGdInLbTObNGReVkT2qq6vZ\nu3cve/bsibQ76OnpobGxMSHrSzh9eMeOHSl/oXZ0dLBt2zbgXGxK+PU999wTceGERVR9fT379++P\njA8LilyitbWVtrY2GhsbI/2eysrKKC8vx+/3c/r0aerr6yeIloaGBq677jp27txJY2NjREjU1dVN\neY87duxg165daQ0ADscHbdu2LRLf1NPTQ2tra05YywzDyH9E5FqgDliK01prD04JmHBpmHLAD1SI\nyH1AN9CgqocSmt/Jkpv7VFVV6Vyqg/Plh88A8Nn3OxK+tLSUgYGB2VySYRhG3iEi7apaNf3Iqakq\nED0wC7kOcoK0rD/biMhngdOq+tVMnTdvXFR9fX3U1tbS0tIy20tJKy0tLdTW1jI2NjYn788wDCMT\nhD87sRiWrCMi24CmZMUNgKp+GdgrIvdOex2z4OQnZsExDMNIHbPgzF0sBscwDMMwUkWwb9Qcwx6H\nYRiGYaSKCZyEEZE/ASpSnKZTVf96qgH2OAzDMAzDyCbXA9tSnKNpugEmcAzDMAwjz2gZhJYhID+D\npPerakr1e0Rk2qBaEziGYRiGkQ6y+I26ebGz7e7Kv0J/biZUxueYN2nihmEYhmHkFiJyjYhck4m5\nTeAYhmEYRqoI4J2FLY8QkZ1Rvy8RkVeAdqBdRLpFZG06r2cuKsMwDMOY44hINVCjqvVxjkX36ClT\n1V1JzLtXVbcmOHypiHxJVT8P7ACagU732AacwOE7Er32dMwbgROuZLx582Y2b94828tJGy0tLbS0\ntEQqGc+1+zMMw8gE4c9O0hWkm6Np4mFhA1QCXXGON0SLHhGpFJHtiYgcEWnA6RWVKL04vac+D3TH\nxtGIyANJzDUtVsk4T7FKxoZhGKmTtkrGJaIH1qdjRckhv0mskrErRspUtS5mf6uq1sTsm9YqIyJ+\nHLFSraobE1qryIPAXarqFZF7w60aRGQdToPNHaq6I5G5EsFicAzDMAwjVcIWnGxv6Vi6yF4RKXN/\nr8bp6j0d1TgdwJOhB3hIRP5/oEpE7nT3d+HE4nROeuYMyEGDmmEYhmEYWaIOR1wcdIOAu1S1eaoT\nXBH0IJCU5UtV75vk0EagV1UPJjPfdJgFxzAMwzDmKaraBdQDbUADcE8Cp5Wpam8a1/BkusUNmMAx\nDMMwjPSQhbTwptNQ9cq5DVieypJFpBFoc2NuaoBqEdk7xfgt01l4krz+bemaKxZzURlGvqEKr/wE\nXmyFjR+ENVfN9ooMw8hSFlXtameLXPYAp2Y6l4hU4riGugBUtU1E1uO4qyZYadzA4rRZblxqgB+l\neU7ABE7eEAgqp86EONEX5ERfiFNnQixfbAa4ecXoIBz4BrTthP5TEBiGn/ydI3CqPw9XbQZPnlX+\nMgxjNiknJrBXVXvdbKd4VAJ+VxgBXAeUich2oDkslJJEZnBOQpjAySHGYkSM89P5vftsiOiM/pJC\n4dLzC2ZvsUb2+T+3wfFfw2hUOYBQAI7sh298BG74PbjrH2dvfYYxn8nROjjTcAAn/ibSmdvNpuoN\nW29cMVOlqk2xrim3QKA/mcKAcchYrZqEHoeILMZJCbseWA+UuYd6cdK7WlU1Iyamuc6zh0dpfXqY\nE30hevonipiVSzz4V/m48RIPKxZ7WVXmYeUSLwuLBZGMCV8jFznz+nhxE83oAPQdz+56DMPIeVyB\nUg1scV9vx4m56XCtNXVujZxunO/0sphqx9U4gcdNMfPWAltxLDrbgaZ0Bh6ngykFjohci5NCthQn\n330PzhvQ4w4px6liWCEi9+G8QQ2qeihTC54puVrJuPtsiKOngvQPjxexRQVw4XIvK5d4WbnEETWr\nlnhYscRLUcE5YWOVjA3DMJIn7ZWMcxRV7QA6gLhWlqgsqsnO3xXvXFVtIkb05BqTVjIWkc8Cp8OV\nBhOecIbnZZpcr2Q8OBKKckuN/3l2aPwzKisVVi7x8rYrirjhkiLAKhnPC3ZeASdfglBw4jFvIVT+\nNnz469lfl2HkMWmrZLxQ9EBGemJPjTyeWCXjXEVEdqazenE0cS04IrINx9zUl+yEqvplt0vovbkm\ncnKZBUUe1q30sG7lxEcyOBLi5JkQJ3pd0XMmyFMHx/jv50ciAseYB9z7XSfAuONbgAfGBqFoIYgX\n3vpJuOXTs71Cw5i/5GcMTi7QM/2QmRH3cajq7lQmdYWRiZs0saDIw9oVHtauOPe4wr2ojHnEig3w\nwX+G9/0N/PJf4Nffhd/aBtdsBV/hbK/OMAwjaWIbbqYT05uGkW8sKINb/9jZDMMw5gBueIviBEA/\n5XYWX4+TyDSj2N7JXFR/AlSksFaATlX96xTnMAzDMIzcJ8suqpaT0OKU+JsrQdK9Ye+RK3aWquod\n7us7gUPJTjjZ47ge2DbDRYbJ6ehqw8hHFGWQx+nmKwxxgEW8k3LqKOaK2V6aYRhZZPMKZ9t9jKRj\nZXOU6CKB1Th9scLM6B4nEzj7ZxJgHI2I5G7KkmHkIWf4ISf4AkF6UQbdfd/jLI9SSAWr+RIlvHmW\nV2kY8xQLMk6V6HThapwaO/GOJUzcWv/pCPrJZOCQYcxHTvCXBDgWETcOQZRhRniOHh6YtbUZhmGk\nSIWIrBOR+4F9qnoGIu6psqlPjU9KetMtBOjHUVfCuQAhS/HJIK/3BukdCFFWar2o5hdT/xGjmat4\nbhhGIlgruBmjqrvd2JtuVb1dRJYAtcAymFlD0RkLHLfFeZeqPhSz/y4Rac01kZOrlYyT4ZXjYzz6\n1DBPHxzD64Vbryq2SsaGYRgzYL5UMs41ROQanIDiQ1Gvu1T1TLTnxw2TSckTNGkl4wQWeVesuIk6\ndqeqfjuVhaWbXK9kPBWBoPL33zvLi68FKC0S3n5VEbddXcziBecsOFbJeO5znHrO8DBKEAhE9gtO\nsccVfJ5yPjZLqzOM/CRtlYzLRA/cko4VJYd8Nz8qGbsWmQM4Xh9wUr8/LyKbcAKKr1HVtEYxzcjH\n4SqutqjXO10/WWRXqgszziECBV7nLQ0pBIIQCM3yooyss4YG1vOfLGErQjHCAjwsppz7qOBxEzeG\nYeQyDcBB4G5gB3C3iHxGVfcB7WRAN8xULS1T1aeiXm/FCQIKW21mHAwgItVATUw3U0RkC47ya3I7\noNYCB9xGYnMar0f4o/cs4uAbAVqfHuaxp4dpe2aYD9xQwh3Xlsz28owsUoifNexkJZ9nmGdZwHUI\nBbO9LMMwjOkoV9Xbo17vEpEHXONIRgIIZxqlesCNwQFAVTeo6ieijietxESk2m3ZXk/8iOlyHAV4\nWkQU6JkP4iaa9at81N6+kL/68BKWLfLwzOGx2V6SMUt4WUQpN5m4MYxcIZwmnqWt5STUPgPkTwxR\nl4gsjtEO9+EUFc6Ii21GAscN/lka75iIrGN8wZ5E52xzrTZTiZYKoEJVRVWbk73GXGHFYq9lUBkz\nZpAXOcXDBOmf7aWkhX6OcZDHGGFiXkOIEU7xQ87wZFazzEY4w0Ee4yzHsnZNY5YRnCyqLG2bL4Cm\n64AZFsGbBfbgfL+3umEuQKSkTBPZdFGJyOKpMqFU9SERuQs4rao/cgOIqsOv071Q95pJCyfDMEAJ\ncJr/4jhNDHEQQTjEX7KczazmYxSzbraXmBRKiDd4mhd4kG5eQvDQQSPn81tcxp0sYCHH+RZv0Aw4\nAWsFLOd8Ps5y3oHHDcxON70c5EUe4lV+geBBUcq5mMvZwiquRWZsNDeM/EZVnxSRjUBVTIgLqtok\nIvvTfc2pYnB2uNukuCJniRsF3T1ZVlW6cONuenDcVWWquiuT1zOMucJzbGWEY4TcIoFhW8ZJHuYU\nLazlz1jB+2dvgUnyU/4/uvkNAYbH7T/KzzjFj1nOG67AOOfGHeEoh9jFEf6BSn6Ih/R2YH+eB3mB\nBwkRQDmXBXCK5/gFnZRzCbfwl2m9ppFDWCXjaXG9P/smOfZkuq831Z8T1SIybd13Ve1T1X2xiiwD\ntKlqk6o2q2oTsExEtmf4moYxJxjh1Yi4GU8AZYQRDmd9TalwhlcniBuHEB6GUBgnbs4dHSRIPyFG\nMrCmIwQZHSduwgQY5gyvpv2ahpGviMjiTF9jKr25AegQEXBSwltx25jHDhSRdTNpZZ4McdxT+4Hd\nQEJWnJMnT1JVdS6Oqba2ltra2vQtMMscP3acY8ePUfWXnwVgZGRkTt2fYRhGJmhqaqKpaVwv6OWz\ntZZ5TgPwiWlHpcBUAqcSJ6amB6gB7sNJ61IcwbMXaFXVw0Ad07izUkFE/EAnToBxWOj0kkR/ihUr\nVpCvhf7isea8Naw5bw3/8UXnnkpLS+fU/RnpZipjrWea4+lnkOc5yf+lhA0s4068LErqfM8U61UE\nmSKg2LGwpDeecZAjjPE0JQwyTDEaZ30eq+OfE8T+8SciM2oDEBdzUSVDnYh0MonhJB1M+jhU9SCw\n2w0kfkBV73MDiWtwhM/ngCYROe2ekjGBgyOy6mOsOGXMIFvLMOYjG/hbXuUfGOYgIUZxXDklKMoy\n3s0qPpjxNSgBemnjDZoY5hDKKB4KOcY/sJR3s4rfp5j1Cc11A5/hGb5ODy+jBFGCeCl072cTy1F6\n+E8cOTMICB6KKWA5F3IfPham4X6UHp7gMF/nLC8QIkgpygIGGaOYAYoJUYTgpZwNXM1HU76mkcOE\ns6iMRNnqxvFeKyLbcOJ409oBYVq96S5gvVuMp81Nz26GSOnlGmBnOhcVZw29rqssmjqcmjmGYUzD\nEm5mCTczyIsc518Y4FlW8kFW8AG8afiyT4QeHuEof4VGxc6E3N97eITTtPBmOhLKNFrGZdzKTs5y\njJd4mNfpYD01VPAOityyICH+hFM8ynG+TiFrOJ+Ps4hrkDRZb3r4Jb/m8xPieQQoZIQCRljEnVzG\nb7OI89NyTcOYK4STktzg4ifBaQGFU4LmQDqsOgkZ1FxrzkER2SQiGk4DdyOim0Uk5SI9IhJ2iW1x\nX2/HEVThujhN7r5enHo4e+dzLRzDmAkLuJQKGmbl2kEGIU4ArkMARUi2oOkizmMjfxD3mIciVvJe\nVvLepOZMlCCDyKR/siteirmaD1PM6oxc3zDmGmHRky6rTlIeQ1Xd56aF34UTfxOuk9M40wVEzd2B\nUwQobtCwqvZOdswwDMMwZpUsp4m3HIaWI0D+VDJOmLBVxzWqvIxj7Eg6IDnpx+FabR6Ktua4Fh4j\ngwyPKke7Axw5GeToqQCvdge5YJk5fHMWDcGrj8LTO2HsDLypHtbdBd701l7JJ7yUMnkws891TSXn\nPurjGM/RwjGeYQNv51JupyRLn/c+St3O7vHwoAQzVlDQMDavdbbdL+ZNJeOEcNPH78YJQfHjZEs3\nTXnSJKSiNzuBShHZgxMAfCiFuYwozg6FOHIyyJFTAY6eCnL4ZICTfaGI8X5RibB+pY9brrQPz5xD\nFX7zFXj6SzB6BgJuO4TH6+DnfwBX/CFc82fzUuiU8z68LOYNGhmiE2UUoQhQynkfq/jdhCv9nuAF\n2vkW3XQRIoQS5Nc8wq95hAuopJIPsTjDrqGlXM/VNHCYr3OG51BCrstKWc7bWcuHKYzf0caYi1ih\nv5Rw2zd8HrgL15ujqrtTmXPKVg046skPXBf1u59z6dl9OJlM9WQ4n32+8MW9fRw5ee6vQq8Hrrig\ngN+6tIgLl3u5aLmPslIhTtC1kQuc7YRffgZCMYXkxs46P5/9G1h5I1z4ruyvbZYRvJRRTRnVDPES\nJ/kWxVzMMt7vWncS56f8E/2cHLcv6Bb2O8yvELzcwh+lbe3xEIRyrqec6xniVV6lmQKWcj7vp2Du\neQ2MPEdEqoEat+dj7LHoomlTdgkQkTIgPP46oCvenAms509wLDO1OGVo1uNYa6rSVdV4Kr3ZixPx\n1wUcdH8+6P7sMLdUZrh2fSGFvjGOngowMgbBEDz/6hh9gyFOnfHSfTbERcu9XLDcR3GBiZycQ0OO\ndSZW4ITxFDhj5jklXMJF/MWMz49XLTj66OSuo8xQwgVczKezek0jB8lBC05Y2ODUtptQWkVEGqIF\niohUisj2KUTOjpjx7bFzJMgunGJ/TwINqVpr4jHV4+hyL97DHBA0fX191NbWsnnzZjZv3jzby5mU\n91SV8J6qEkKqnOgLcfRkgCOnHHfV04fGePyFUcCxhtZcU0xx9z5aWloYGxvLi/szJkcJoQzgSbLo\nnWEYydPS0kJLSwvMwSDdaFS1DWgTkQbiF8etjBnfISJx69q51pvumN2NOFohWYHTAWzLRA+qMKIa\nPy1TRO5X1c+5v2+CSAWuHpyI5jNRY+9Md4GedFNVVaX5XulXVTk9oBw9FeDBxwcpK/Xw2fc77TxK\nS0sZGBiY5RUaDHfDgxWgAQjEPA9vifPzvb+E8qsju0OcYYBv0s8/EaKHYu5gEZ+icPznjhFFKzs5\nwQsT+lEJHjz4uIr3cY1TccIwpkRE2lU15VInVatED2S+XuYE5B9IaP1hgaOqdTH7W3E8NtvcmnPV\n7rgJZVjcci7twMZwCRcR2YJTtiUpl4KIbMuE1SaaqSw4keJ9qhrp/hku7ici5TgurNM4VY1zWuDM\nBUSE8oVC+cJCHnsqXqNBY9YpXgYffA1e+Xcng2qkB0TAUwhX/TFcVueMcenlzxngX93mAkMADPMD\nRtiHl4vrz/kBAAAgAElEQVQo558p4LLZupucpZrPcZxneZZHOMmLbufwEBdxA1fxXspZO9tLNIx8\noQ5HtBwUkZ04MTVxa8y51p2aqPp04Li/2pK9aDxxIyKfxYnJaVPVT4jIeuDamRpQpmrVEDf1LJwm\nHrWg9WCpAoYRoaAULr/PETPHf+wEGF/4bvBM/O82wAM4MSPRhFCGCPAKwzxmAicOgnAeb+I83sRZ\n3uB1nudCqig2954xW+RpFpWqdolIPY5QacDpVDBpEV3X5QVEXFZ3A5tSXYeI3A+cAm7HdZtFFRme\nkZco5cehqgdFxCoKZ4nBkRC/fHmUYz1Bzi+3Ojg5jQicd1sqE6RtKXOZRaxiEatmexmGkRWannG2\nKFLqhi4ijThBvk2ue2qviOxV1a0JnL4bp6dUx7Qjp6czbNVxDSfjljmTCeMKHBFZl0xdm3CsTipz\nGFNz6ESA/35uhF++PMJoAC5a4eUdlcWzvax5T4gxzrKfUq7El/ZYxRCjnJ5+WJKM0s0QnSxm4xSt\nBgzDSIosWXBqK50tctkGZtwN3Y2p6Q03slbVNldcHBSRMreDwGTnNgCN0RadFIkOXo4VNDPyEk32\nOBqAe2YyYZrnMID/7Bji20848RlvWlvA5utKWLcyD22hc4gxTnGCPbzBNyPdrMu5ndX8Hgu4OOF5\nink3w7Ti9Ggai+wP4SUEvMAP8fEq6/koK7gpJUHSz3Mc4185zU8RfHgo4jx+h5XchY/FM57XMIy8\npRynaG8EN9D4walOcgOLW8PiRkQq02DFqYgyjES89m4BwLQKnKVusNFMESwuJ21ccWEBv3l1jN+8\nGuC5o2MU+OCWK4u59HwfHiv4l3VO8BBH+BIgaFQn6W5+QA+PUcatbODLCc21jK8R4DD9PMAA3yDE\nMCOUcIqlnGURzn+lF/g1X8TLAn6Lr1M4A0vR89xHP08TYhQnxmeUEIO8ShNHaWQDX2IZqbjTDMPI\nQ4PoAZz07kgrBDeupjdsvXGtPFWq2uS+rsYRRm3uWHAaZacqcJpw3GNL3et04TTW7lbVO2Yy4WQC\nJxHfm5El1q7w8cfvXczrvUF++twIj78wQnvnWW6+rJDfvW3hbC9v3tFPO8ponCOOJaef5MoR+FhL\nGTsJ8G5e4M8YilOkLsgQIIxwckYCxxE3EzPvwvsGeM4EjmHMQVyBUg1O3QQR2Y6TpdThWmvqXHdT\nN066eFlM0b5qHG9MkytoWt390U22U47DdROYbnfL0oSdcPenUicnrsCZLIPKmF1Wl3nZevMC3n9D\nCfd/+wwnz1hF3LmEUMwYJUD/pCMMwzCSwXUddeBUDo53PNxuabLzd4XPda06Gf0gUtV9ItKZjhje\neRPIka1KxqpKMASjAWU04P4cU0bCv0e2c8dGxx2L+n2qYwG47HxfpBpnrlYyHnnjDY5//ess2LCB\n5e99Lx5f1D85DULvD2DwOVjxUSg8b/zJo4fgzDeg5AZYUO1kJWWBUHA/oeAP8PjuxOO5evoTIiiF\njFLCGwQ4iI/YRIAErh0VhxNLkGH6eI6FVCBRnzFBAhzmCYbpo4JbKGK8VW+QLoJRViGTSYaRgUrG\neZomPlu4/S7Lo3b1RBUQrhORBzj3cRV9LPFrTFbJeK4RrmQ8MhYlMsYcwTASV1yc+30k5nVc0TF2\n7nVoBm+pzwOFBUKhDwp94m7RvwuFBeOPXXVRIRWrnf9RuVbJuG//fg41NND9/e+jgKegAE9RERd9\n5jNc8PF7KAh8G45/GYIDoGOAwJLbYU09eIeh+34Y+qnTt0kKwVsG5Z+Dso+BJ/1uOdVRQsGHCI79\nFU7phTHAh3iuwOv7UzzezYg47/UgL3KILzLICygBhAAlDLKQQbyEXPFRQCFVLOLTFHHLOEEyGUFG\neIG/43VaASEUie8RnJg7D16K8LGI9XyEMm7kFf6b3/BDt81DCFDWciNX8m6UI7zKvzLIy65wCuFU\nkAkiKB4KAGER17Cez1NixfGMeUjaKhmfL3rgvnSsKDnkfyVWyTjXEJG7cJKRluL0uWwN17oRkcdw\nXGPgFBOuV9WvJn2N+SRwmvb8nAce7SfVO15QJCwuERYv8LCwOL4AKYoVKAXjxUqBD3eM87vXk9rf\n1bkkcN7Yu5fnPvYxQiMjEBrvRvOUlPCW7wxTUF6ChAZjzhQo8kKxB+LFuEgp+FZDxStpX/PY8LsI\nhf4biPceluLx3klB0b+N2ztEF6/zNYTd+AghcWJnhBIW8kkWJ9GmZYyzHOP7vMK/EHKrG0+kmOOU\n46Eg0kX73DU9rOAkSxgaFwQdPUKANWxhDR+mmAsSXpthzDVM4MwebuXiptiwGLfo305V7XOFUKeq\nPpXs/PPKoLZhjY8P3FjC4IgyPKYMj07+c2SMSYXQ4IgyOKK83htCBIoLxNkKheLC2NfO7yXuzyJ3\nf/h1eFxJEXOmO/joG284wiY0MUYoNDREwSLiiBtw3vHA5BPrAARPpG2d46bWV4kvbgAGUD06YW8J\nftbzBY6zh9CE/nPuvAwR5HhSaylgEWv5bYY5yRHiZ2sGGQV0grhxrhnCR2ASceOMKOJC1ifdG88w\nDCM9iMi9qjpZuml3WPSo6kMici9gAmcqFi/w8M7KkoTGhtRxS4UFz5ArfkbCr6ME0UjM6+FRpW9w\n/OtE3FZLFgiry7ysWepl9dJzP5eWCmLp4IZhGLlLlmNwWn4DLS8A+dsNfaovtaaY1xUzucC8EjjJ\n4BHXGlMoUJraXKrKWJC4lqIhVyANjChv9AY5fjrIL18eZWj0nCIqKmCC8FlT5mVVmQdPiq6tTOAp\nKpoyIDgUAI/XaS85AfHg/Luf6O4BD0hRupYZc90FoD7iW5B8CAsm7FUdhdG9UNDnLDnuLRcgzKzi\ntIdCBB8aZ03hqJzJCCE4RTnivY9ePGTofTQMIytsvtzZdu8nX7OeJxVmcTK5y+IOnAYTOFlA5Fw8\nTiL1YlWVM0PK66cdwRP++cJrYzzx0rnYlLdeUcRH356i+soAqz/yEYIDAxzatYvg2bME+/uRggLE\n62XJzTczIHexeEkLnPkxEAIdBc8i8BTD8k9CwSj0fgU0AHoWxBUIpe+E5X+ekTUXFH6LwNj9hILf\nwJEPg4SVrcf3cXy+7ZGxGuqG4b+HkX8EDbLMO8TZ4oUMF3oBL0gIoRTwspB7WcgnZrSmdXwID4Uc\n4UGUIEGG8FAIwEpuxM9NvMwT9HCQEEHXNVVMEQu5iHso4iBv8DAAIYbwUAwoy6nhAn5/xu+VYRiT\nkH+F/maTZUmMLZ9+yERSCjIWkduAAzNJ38o24SyqfCOkysvHAux/ZZSOrlHODilFBXDNukLuuLaY\nC5fnZhYVgIZCdP/whxz+27+l9IoruOjTn2aB339uwMgReP1/w0A7rP40LH0PiPsJoWNw9jtw+p9g\nwS2w9H84AcaZXrP2EQz8C6Hgg3i8H8Pr+x1ExotI7a+D0X+FmPiXoAgDRSWMFJay0PcVSngv4gqS\nVAgR4AQ/4SjfZinXcCF3UhT12dDHazzH9xikhyvZzGqujGRtBRnmFD/kdR5hObexivdZWwbDiCJt\nQcYXiB74VDpWlBxSn59Bxm4g8QPT1btxdcbt8XpeTnuNFAXOg8D26AWKyOJcFDz5KnAefXKI5l84\nmTSLSoS7b1pAZUUhhb7x/pBcFDhzFe3/KIz+++QDZDWyNLnAYsMwZoe0CZwLRQ/8z3SsKDnkM3kr\ncJYAbcAWVT08yZhrcVLIN85EV6TqonoAKHML9oSpBf46xXkNl+suLqKnP0R75yh9g8q//2SApw+N\ncd2GQq5eW0CBL/dicAzDMAxjKtwU8DrgSRHZg9MCoss97MdpD1ENVM/UaJKqwGkG9nMu5lGA9eSg\nwMlWJeN0U77Qwz1vWcDbriyi9alhfv7CKAc6ne3mywpZNvDjnK5kPGP6u+Dlv4NDX4Ol18KlO2D1\nHW4Q8uwQoJceHuRk6VN4S9azYqibstGzeMaF+5aCZ9WsrdEwjMSwSsazj6p2iEgVcD+OngjrCNzX\nVepUXp0RqbqoNqnqvph916bSHCtT5JOLanAkxMvHAxx8w91OBCNZVSWFwrqVXtav9HHjpUWsWerE\nrMwZF9XZl6DjPuj+hdPKQd04F99C8C2Gq3fC2o9mdUlB+nmNL9DHD3E6iDsNKj0KEGLZ8BlWDQ0i\nnuVQvAOKPoLIxKwrwzByj7S5qC4SPfAn6VhRcsgf5Z+LSkTWxYu9EZH1AImImsnmiCZpvRkTY9Mp\nIutihtwN5JzAySfu//ZZjp8+l967qER4T1Ux120oYvVSD565XBOnazec/PHE/YF+Z3vyk1kXOP08\nQR+PTiicFxIAD6eKyyn3/ROFvi1Wr8gw5jOWRZUoDTguqHEkaa2JO0c0MzGoHRKR29yyyW1AO+Mr\ngFwL7JjBvIbLR9++gCe7xug6EeDIyQBnh5TvHRjmJ8+NsH6lj/WrfPhX+diwxjch2Dj/yc3WIU49\nmkmOSTFScGNC/aYMwzAMlorIzhTOF5weVlOStMBR1eh89I+r6k/GXdWJejZSYMOaAjasKQAgEFSO\n9QT51Suj/PT5EZ45PMYzhx23zU2XFfJ7t6W/8WRGUYWBF6H0klmNp5nASBcUrAFPYpWux5Obosww\njCyS7UrGzzgb+VnJeGs2LpLq4/gDYJzAycX4m3zk8IkALx0P8MrxMTpfD9A36HyJFheAf7WPDasL\nuPGS1GusZI1AP7z2dehqgNETUFAO6+vhgt+DgqgkvPLrnIJ/IhCMbjTpAW8xLLspfWvSAPR9F058\nCYZ/DVIAy+6D5Z+Cwgsjw4pYD4CwAGV8Dy0PC/CyBO/0f0wYhmGkjc1vcrbdP8u/SsZxKhVnhFQF\nzgG3CM9SZtjt05jIi6+N8dePnAVg+WIPl51fwIY1Pjas9nFeuTcn2zNMyZEH4DefcUWLGwg9chxe\n/jy8tAMu+UtY/8fO/gvvgZWboKvRyaIKjTi9Hda8By6th/I0xdIN/BwOvg90BELOe42OwKn/Daf+\nEZbcBWu/CUAxFVzBz+nl+5zgAQKcQAlRwuWs5BMs4m0IOWSNMgwj+1gWVc6R0uOI7gQqIkvcwn+d\nqmoxOCmwdoWPZYucL8y/uGcJRfneZfzYtyBe9/Cgu++1fz8ncACKlsPlf+oImpP/BUveBMUr07um\n/h9DsJsJ7iV1W2H0PhgROAAeiijnTpbyAYZ4Bg8LKZ5Z/zfDMAwjC6T0Z6eI3Ckit7nCpg2nJk5s\nF1AjSYoLhY/dWkr32RC/eHFk+hPmKh4frKpOv7hJAUFYwJtN3BiGYeQ4k1pwEmy58FVgD7DTYm/S\ny9oVTr7hWLxm0POF0FkYegQKq6DgsrRNq4ToKTyFLipj2dnTac19OstxTvISF1BFYapt6A3DyB8E\nSxPPMaZyUe1g+nTvbar6UBrXkzHyrZKxR5yk41eOj7HpTUWT1r4JV+PM6UrGa/8Q+n/jxNMEz57b\n710IngJYF9PAZewVOPs3MPhvgMcJBi6shMWfh+J3zjj7KsAgx/hPDvFNAmVnYPHleENB1p04ynk9\nb1AQCoIsABSW/UHC8yohjvMMz/EQ3byC4OVXNLKOt3A572MJ589ovYZhZA6rZDz3mbSSsYjsB+5V\n1aezu6TMkE+VjMN8/8AQD/9qiFuvKuK337JgyuDinK9kHArAiRbo/Cs48yQsugoq/hRWfcAROWH6\n/xlOfxIIEtutG1kIvstg9f6kLz/MSX7Oh1GUkFuJOIwnBGiQ3+p8lQVL/xjKfx+8iXfcbuV/0UMX\ngZh5BQ8evFzDR7iMdye9ZsMwMk/aKhn7RQ98IR0rSg75nfyrZJwtptKbG4AOtzJrG04jrLZ4mVKJ\nlEw2kuddG4sZHFEee3qYjq5Rbri4iBsuLeTCZd78q5jr8cHqDzjb2GkomCStOvAbiBEKEbQfAi/O\n6PIjdCN4CDJRBIY84NPFjFz8CAtmUMapj6MTxA04lp0gIXo5MqM1G4ZhGONJRm9MJXAqcTp59gA1\nwH3ALhFRHMGzF2h125zXkabqxSJSDdSoav004/aqalaKBc0WIsKWm0qoWO3jFy+NsO/ZYR57epjz\ny73ccEkhN19WxOIFeZiePJm4mU1EwCoRG4aRCtks9NcOLR1Afhb6Swi3FdR9nEt3FWATcF0i50/6\nONyeELtF5C7gAVW9T0SW4IidauBzQJOInHZPSUnghIUNjrDqmmZsA0479TmPiFBZUUhlRSH9wyEO\nvDLKEy+O8u0nhnj5WIBPvWfRbC8xzXhxkvtCkxyPL+gGOMpRHqKY1ZzPuylg/PsiCKFYl1cUAYY4\nzFMs4nJ8FCW14qlr4IjVyDEMI+1s3uhsu3+Uf4X+kuBzOMaUaMoSPXlavamqD4nIehG5E8dF1YzT\nxpwowZNKT4nwddqANle8THoDIjIvhE08FhZ7ePtVxbz9qmK+/PAZRgJzsEXAwj+E4HEY3OtYVXQI\n8IAUg+9iWHLun5qidPMrDvJvnOEFlBCCj052s5pq1vJBFrIOgEVcjJ/f5zD/lxCjBBly5/ASQjnD\nYl7lZzzJ41zGLVzFO1jE8oSWfDN/xJN8gz5eJUQAJYSXQhTlQq7nSu5M97tkGEauYUHGmWCvqu6L\n3iEiCQfTJvQ4XGvOQRHZJCKqqj9y9/cBzSKSzQCnapx4oOosXtPIFr4LYNm/Qdnfw8Bu6P8KFN4E\ni7dD4fj4mJP8lGf5AqGoLt9KAIBjPMrrtHIL38PHAgQP6/kwa7mHU/yc5/k/DNBNL0sYZAHOp5NT\n5O959vESP+NjPJDQklfzJt7JLk5zmOd5hBM8xwaquZg7KCbxYGXDMPIYSxPPBOp2S4j26iQcEpOU\n3lTVfW7F4rtw4m/CdXIak5lnprhurAcBixif63jLYXG9s01CgMEp3D9BwEuIUWBBZK8HHyt5G6/R\ny8s8RLxGmSGCjJF8gcWlrOVmPpX0eYZhGJlmqvhWEamNelmmqrummWs7jugoB1DVTBX4bQLaGR8g\neS2ZEDgQsdo8FG3NcS082aBMVXvzLoPIMAzDMGaB6eJbRaQhWvSISKWIbJ9M5IhII47rqM19vV1E\ntrjhK+mmLo6LKuFU11SiHzuBMhHZ40Y6Z5RU38CTJ09SVVUV2Zqa8rejRCCovHbsDdrb2yP3MzIy\nknf3N8bTnOb3Ocn1DPAvhKJSuEOc5QxNHONGuvkUY0xMD/dSgk4RjKwoHgriHi2gCO8k+l7w4iOP\nOrUbhpEwTU1N4z4rIcFgu+kIx+Bke5sGVW1zBUzHJEMqY8Z3MHWWUm1Y3Li0kaYs6ljiiJvbgI2J\nnj9lqwacTCU/zs36o7ZwEHAfjiKsBz6RzMKTwQ0s7k1ljhUrVpBvhf4AgiHl+Okgh04EOXwiwMET\nAV7rDhJgGTffuJr/+KJzT6WlpXlzf8M8Sj9fJEAnMIIjZ/6cfv6MIu4iiDJICyAoQwQ4wiCPUMAV\nlLGDYt4GwEreytX8BQf5N/rpIsQYHvef9Erexjo+jG+SdgmXs4liFvEk36WfbgKMuqJGuZibuZp3\nZeW9MAwju9TW1lJbe84jIyKnZnE5OYGI7MXpTNDrWnz2TDKukonfxV3EiKQ0rutaoAFY7+7qw+l5\n+dVEzp9K//XiBCh0AQfdnw+6Pzuy6JYC583zu28uOIKrzPUDNqvqlGnl+UQwpLR3jnLwjQCHTgQ5\ncirAqBM3S0mhcNEKL5veXMy6FT4uvyA/Q/Z7+TC4wcDnGECBfvYQwsP42JggSpBROjjJvVzIS4CT\nnr2St7CSt9BPF0doppjVXMD7KJymNIQXHxu4iQ3cxAk6eY42VrKei3krhZSk8W4Nw5g35OdHch1O\nnMtBEdkJdE3hLfHj1MabgIiUqWpKhog4bFLV20XkLjejewlJxOBO9Ti6cJRTD9kXNOOIfbPdgCj/\ndIFQ+cjzR8fY3Tq+2u555V7efmUR119cSGnxXKipMplbCRxhM1X6e/zuowvxcwXbZ7SalVSw0rqD\nG4YxD1HVLhGpx4nTacApAzOTcJByUvS0xGFcE29V7XOLDSfEVAKnWVV3A7gBxZvc/T049XAincZF\n5E5V/XYSi56Aa52pBra4r7e71+mIGVcLbMWx6GwHmjKgGmeNqy4qoP4Di+hyLTiHTwY41hPkWz8d\n5Fs/HWTlEg9rV/hYt9LLtf5CViy2vETDMIxZJ0tp4k0/cLYoUoohcoOGG1S1yXVP7Z1hp4C4lp1U\nEZGXgdtF5DHgMRwPzo8SOXcqgROpqBYd6BMu7ici5Th/ap/GqTaYksBxhUwHMKVVxk1Hy/0I2hki\nImxYU8CGNeeCYweGQxw+GeTQiQCHTgTofD3A/ldGef5ogE9vzr9Kxh4uQOlB6Y85UuCmfRfgWHnG\nW2uEUnyszdIqDcMwco/adzlbGLmDGccQhWNqwmEeqtomIutx3FXxXE6R1PAowqniaTc0uNrjYnet\ndZzropAQU7VqiFv+OZwmHn7tvhk52Fxo7lBa7OGKCz1cceE50fPlh88wFszPSsYrOMAw36WfvybE\nERQBQpTwIUr5HyhjnOErDPCwK3gUHxezhE9Twh2zvXzDMIyJ5Gcl43KcjOgIbqDxg/EGq2qHiMR2\nGihj8gytlHA9R/VAu6ruEJEukui9lfLjUNWDIpKJ/HdjjiIUUsIWStjCGB2M8RuKeR8eFkbGLOMf\nWMoXGOR7FHINhVw5iys2DMOYkxzAERARr4grYHrDFhnXylMVVcyvKaZsyz1krtiv3w0y3gSRYsO3\nJXpyWvSmqiZsMjKMaAqopGCSDEMPS1jIh7O8IsMwjLnDVPGtrrWmzu0B2Y0TJFwWU+24GkfENAGo\nap2INLgFd8uB7gxWMg5bl6LdFX5SicERkXWqeiiVVaVjDsMwDMPIC3LURTVdfKsbfzNpTxw3W3lX\nzL7Je+iklwoRqYGIZakGJ6U9ISbLOd6YjBkoFrdXVU51/e7r66O2tpaWlpbZXkpaaWlpoba2lrGx\nsTl5f4ZhGJkg/NlJEjEd0+KdhW0O42Zyt+HE+d6OkzWdUJE/AFGNH6jqVhCsAx4Mdw+fdjJH2FQD\njar6VKKLyAZVVVWaL5V+E+HLDztZ+p99v9OturS0lIGBgalOMQzDMGIQkXZVTbmBc9Vlogf+OR0r\nSg55C2lZf76QTFmaqbKongTuE5G73IjqJTgpYp2cK+ZTBlS422lgj6pmrGWDYRiGYRjzAxH5LE78\nT9gSIzhtG1ITOGFU9SHctHDXquPnXB78QWCfK4YMwzAMY36S5Riclp9Cy8+AdLrYco+uWOuU6ylK\niKQehytkTMzMdUJBePIHcKIT3vI7sGjZbK/IMAzDiGLzW51t9yPErVk3R4hXPDDhIOMcjPk2Zo2B\nXvjxV+F7X4bRQUfo/McOuO5O2Lwd1r55tldoGIaRm+RoFlWe0yMie3DCY8B5lzfhtGuYFnsceUQg\nqJweCNF9NsTZoRCLStLYeDMYgD9aB4ExR9xE88Qe2P8d+OS34Lr3p++ahmEYcwUTOJmgjomtmWIr\nKU/KjB6HiCwG7sYJLi7DMSN14hQPOjSTOQ0YCyo9Zx0B03026P4McepsiJ6zIU4PhIhOelu7Io3/\nm4IBGO53rDaxhILgVeg9lr7rGYZhGMbU7I3uhQkgIgmnQyf9DekG+ISL7XRFHVoKfE5EHku1s/h8\noOuNAE91jdLdH6L7TIju/iC9A+NT9j0CSxd6WLbIw6Xn+1i2yMOyRV6WLfKwfJGz3zAMwzDmKOrW\n5IvWGnXAjkROnpEJQFXvm+xYMhHO85l9zwzzq5dH4x4rKxUuOa+ADat9rFjiYfkiL+WLPBT6JHML\nkinmPg+4aRTKH4TRGii8OHPrMAzDyFfmeOG9WaAJx5gS/QV1LRkUONO1sM7JFtfhSsabN29m8+bN\ns70cPr6plPffUEKP64KKdUu1d45OEECLS+ScFWexh+sqCnnml/9JS0tLpJLxjO+voAg+9yg89BfQ\ndQBCY3B5EN4msETBpyA/h4NvguLrYPmfQ2lNet4MwzCMLNPS0hKu/D6X06xzFhHZqarTCZW6OC6q\naxO+xmSVjKdY1CYcE1EnTnOuMMtwauTsyUUXVb5VMg6FlN6BEN39IU72OYLn2SNj42Jwbr2qiA+9\nrRRIcyXj4y/B45+Gyx4FXyj+GCkEfycUXJCeaxqGYcwCaatkfJXogebpx6UbuTw/KxmLSA/wAI5m\neDoT10jaguO2Kz+A05Ihut/UAeB+VZ3LOflZw+MRzg4r//XrEZ4/Okb/sCKAf5WXqy4q5KqLCli7\nMkP20DWXwB11cPxxCJ2ZZFAh6Ehmrm8YhmHMB+7Gid1VnJ5TjwEPpStZaaYxOH241Y2NzNH+yjk3\nlQjcfFkRN15SSMVqHz5vBuNxDMMwjOTIdiXjfdDyYyB/XWxNqvo5t0v43ThGkz8FdolIL47g+Q9V\n/c5ML5C0i2raCZNohJVN8s1FBaCqHOsJ8usjYzx7ZIxXjgcIhqC4AC6/oIDbrylmw5oCIAPNNgcf\nh6PVoGNATOq4ekBD0HcbXPwXsOwtIEKIUY7zX3TxTQIMsp57uIB34sNxo/VxmBf5Nq/xBGuo4lLu\nZCkV6VuzYRhGkqTNRXW16IEZfxXPHLk4b11Ud7mtoGL3+4EtOIKnGnhMVd8xo2vMIAbnfiZXjAJs\nVNWEqgxmk3wUOACn+0N0vRGg640Azx8d49Xuc2Jj09VF/PZbMxCDE2b4KejeBf3fAVUIBYAg9Ar0\nKIwJeBegxat4+cYtHFn8axQlyBAAXopRlKVcTy9nOMurhAighAAPXnws5DyupZYVXJXetRuGYSRA\nWgXOd9OxouQQf34KnERxNYcmEJA8gZkY1BqBemDvDM41EuCZQ6M8/sIIXW8EIrVxfB64aIWX6jcV\n4V/lw7/aR/nCDNfBKb4Gzv8WBE7Bc++Annbog3ORzgrBAQY5zqHSJwgxfj1BhgE4SgdBYtcaIsgo\nfWeic2kAACAASURBVBziKb5KDX+f2XsxDMMw8g7XjXX/TM6dSZDxQRFpjU3dCiNT1VMxEuKXL4/S\n0TUWeb2gSPitSwq55LwC1q/ysTTTwiYW33IIXA29k/c4kyktgdP9m8jJygKGYRjGLCIi69yA4/0z\nOX+mQcaTBhhPJnyMxLm3upQP3FBCp+uWevbwGPueHWHfs07W0tJSD/7VXm69qphLzy9IfOLh00AI\niuN0B1eFvi5Y4o9b9C/gG0V8gjdgYsQwDGMC1osqE3SJSCNJdBCPJuHHISKPquodUxy/H6gE2mfi\nK5vrBIJK/7AyMBxyfyr9IyHnZ+z+4RADI87voTh64vRAiI7OECuXeBMTOKeegSd3QaerS9dthsp6\nWLkRRvvh+a/Dr+6HwROw6AK4fgdc9iG0oIQhOujhKwxc/QRceSGLXx2i/OXTFPUH3Mk9FAZ8eEM+\noDgSfxPGSwk+QihFqBuhE0bwIPhYzpUze1MNwzByBRM4meBuoAp4cCYnJ/M4mgBE5BqgNzpPXUTu\nxQkCul1E1ovIvar61ZksKFNko5Jx99kgrU8Nc2bIFSwj5wTLyNjk5xX6oLTIQ2mxsLBYOH+Zj4XF\n4r72UFok7msPC90xJUWCRyRSjTNuJePRM/BINfT8GoKjoK646HoIDn8fvEtgsBfwQMANUO7rgv/6\nn4zt/wxHP7SBseJBlGHwKHiEvrWlnLmwhOLTAS78RQ+elZspuKSetxdewwl+RhffpJ8jCB4KKcPP\nhziPaobo4SUe4RA/QvCgBLmIW7iED7AYKxZoGEZ2sUrGuY+qNgMzLp+YcBaViGzDCS7uxekg3hm2\n6IjIozgllQ+5r+Omf80m2ciiev7oGF/54dkpxUwYrwdWLvGwcomXpQsd4VJa5BkvbMJiplCmjW2K\nm0V1sgO+cwuM9cc/adRN947D2YoSjr1rFVoQ/9+HaAHrR/dQWFQ58Vy6CDBAGVchMfE3AYY4wbOs\n4EoK3PRxwzCM2SJtWVTXiB5oTceKkkNW5l8WVVRsTUbnSMaCU6aqG6Imj7bULIu5UG9SK50jXHFh\nAf+0rZzRwCRupzhuqeOng7xyPMDgiE4aausRKC0W15LjCJ/SYuG6DYVcdVHh1IuSKQKShSnje0Wn\nOCyFUBQnlgdYNK7A9Xh8lHAe109+UcMwDGOus1FE/Kr6o5mc7Db1Pg0cmmpcMgKnM/qFm03VE34Z\nM3ZeR6IW+oTyhZJUGncopAyOTiGMhpWBEWd/99kQzx0N0n0mNL3ACQUnPzaN9U6nNBrFn1dRTvM8\nAQZYQRUyIT3cMAzDSJWWR6HlMSAPXWyq+pCIXCsiDwAPJip0XGFTDTSq6lPTjU9G4Fzvlk8O+3nu\nBsICZ2nMWD8wI2U2X/F4xI2vAZi+x9SXHz4zvYosuxRW3+j0lNKg0yEcQHzgKYCl58HZE4COd2MV\nLKKkp4CiwPmMFJxCia5mXIggLOAmfKyOnBJklGP8mFf4FiN0A4KPYvzczYW8kwIWJvpWGIZh5CWa\nofaA8XjPu5xt97+Tl/0fVfXJ/9fevYdHVd0LH/+umVzIhZCEhKvcAgh4q9y0VlTUwFEr1ktQ26da\n2x5CPee0ansE7fHVt32PVdS2Xtpaora2PXosYL1ErRawaqFaCZeiKCAkXAIIITdyT2ZmvX/sPcNk\nMpnMJHtm9kx+n+eZh8zea/asNXuY+c3aa/0W8B2l1HVKqVUYgVolRmeK9ypQLjDZvNVjLMx5a7jP\nEUmA8wDwFEYK5Xrzfr1S6gFgufnvSozoqq7Xo4jYSc2Cr6wzBg7/8+fw6W8ADad+Hc7+AeRNA1cH\n7F4F/7gfGvZAwRlwzg9JmXINE52ptLOTOspo4nUghTy+Rh7fJJUxvqfpoom/cjMeXN1mUblpYzfP\nsotnmccvGMqk2L8GQggRA1qBW2ZRRcwcr/sigFJqJkYHSb65uwpYbwZDEQv7dJgLbF4fZJe3YluA\nu4AKO65FNagNK4ILn4DzVhiDitP8elNS0uG0m2DG16HtOGQWdnvoEKYzhp/h5kcoUnEwpMfhOznR\nI7jxctNBClm08rkEOEIIESdKqWJggdZ6ecD2FcADWuuwx84qpUr97uZqrR+yoo5mINOvYCaYiONN\nc0r4AvPuJq31I2bFqoDvWFUxEQWpmb3vU6pHcOPPydAoVEgIIUQ0eQMbjDx1lUGKlADLgszUXaq1\nLgtyvBX+QZJSapZSaplVQY6VIkn0Nwxj/aktnByHM0UptQm4VGt9Igr1G7Q8WtPWoWnp0DS1nRx4\n7B10fKzRzYhhEVzwrT0OvysDVxd8YymMHNX3Y8Km0Lh63evBRd/LNQghRAKz6SUqrfU6YJ3ZU5Mb\npMgajOEl/oIGN6ZuuUG01luUUrZM7hvJ6ViitV4YuFEplQuUAo9YVqsk4/ZoWjt0twClqe3krKjm\nNiNzceDfvU1y8k4bDyuL8fat8MRD8MbLRi+N1vDYg1B8Odx2F8wa+JTtTEZzKt+ikhdw0+m3mngG\nABO4kgJ65ssRQggRP+b390qtdaXftlKMMbahHrcaIyZoMHuI/mhxvXIwhsRMxgjKGjAGH6+LJH9O\nJAFOVbCNZgOD7rOTWGQyBli/vZ0dB7u6Te9u7eh9vlOKE3P2lJHfxj+L8VC/ZH+9Jf4Lmcm4cg9c\ndp7Ra+MJSOj351dg3Z/hzffhzLMH1GaFYjKLKeJajvIBe3geFy1M5kbGcDFO0gd0fCGEsJrVmYy1\nApczHmkxgidrDYc57sY39kYpNQuo7GM8zlKMtaGqzMlFlWbGYUuYU8EXmM/hf0ktD7hLKfWXcMf5\nRpLJ+NreDhpqn13EIpMxwO/faeGfVZ20dGjcfbzv8rIc5A89ufyC/1IMJ/8+md3Y6ej9Mk/QTMb/\n3AJXXwxNvVw9zBkGv38Z5s2PrJFCCJEkrMpkPGu20n/7e+wDnOwhnrDq771EpbVeGqLMylD7/cqV\nYgQhJcAarfXiSOrcx7FDroQQyUoJkfTgTA6WGlkpNRGjG0kAN8/PgvlZaK1p69S+Xhz/Hp1m33ga\nY2xNbZOH/TXGNleIvHwZaapbAHTetDTOmSq9I0IIEW9aKdwp0R+E89un3Dz7dLdfzwVWHFcpVUIY\nq3abq3uv0FqXmZenViulVlsY5PTV6xJ2IuFIpok/rJRapZSaxMluoyKgLtQq45EKMZXNO9anATOg\nCixjJ0opMtMVmelQGGYHqNaaThe+AKglIBg6GSR52HOki06XDh3gZGZBZ0fv+zs7jDIB6tnLTl6k\nmUOcytWcwjychDHeRwghRFR9c4mTby45OcFkWErncYsOvRRYEaqAeQmrwTtmR2u9zowJqpRSuZFM\nNQ+h0Uz8txeo9ds+HCPmCHu8T0Thptb6erOBl5qbHuxvAp5AYUxluztgatpmu05N6y+lFOmpkJ7q\nZHgfs7IffjmMSWtTp8ELf4ZHfwIfbDDH4WhwOmHmOXD73TDzZM/mIT5gB8/TzBHcdAEetvAkW1nJ\nZL7MdK4jlRBTzYUQYhBzO2OYyth6xRhBTij59Fy2qcEMSCyhtV6vlKow6+O/sGEFRswRdubmiPvT\ntNZbMKaK+yil/tObD6e/wpjKVqqUqvULaCqBuQN5zkHhgouN24F9sPIxY8Dx0tuhaEq3Ym66eJ8H\n0QED1ly0A7Cbl8mikCIui1XNhRBCxIBSqtcVks1OjTnmtPEKYDlQ5rc/F6NXx7JFts0gJqxxNqEE\nDXCUUndiRE/hUMBsoj9NfDbdl4CIqKtq0Bs/Ee7/eYgCfV3WVHgGMFpfCCGSmUbhDmMdwVgzA5Ri\njAHBKKWWYUy33hKkeLBlloqBG4Ays7dmqdkRUYsxZCQ3lsNFIpnU1FsPznDgIYJfKgom6o0LmKc/\nCwvTQwuDwjO4l4EXQogk43fVpdfvS/P7Neg0XfN79qGAslH7zldKPUjvU/e9HSoDCnD+GMnYGjOa\nizqzK+x6YLF5EwOk8dDGeup5nJlspplhHGEEJ8jB+35XpACaLHpfykEIIQYzjcJlwx6cBLQSI4Ba\nPdADBQ1wIh04bK5DFXXmNb4yoMwcZLwyRDpp0YdmXqOO/8JDK5oWFJBNI5Npxo2TA0yhiTwmcimn\n8hWyGR3vKgshhEhiWusqpdRarfX6YPuDrJnVKxuunBGcUqrI/zIVRpS3Er/BTqHU1NQwZ87JGUOl\npaWUlpaGeIR9Ha5zs/9wA8eP7GXOf98JQEdHR8Tta+BR3NR026YAJ26cuJmCh3H8nhRzyQUhhEh0\nZWVllJV1+9qwJI8MgDtxvlJtLVQiv94Cn2AS4myYU8jXKqUmBwQ5hDv3vrCwkFhkMo4Wj9Z8tL+L\nt7e380m1ixRnNt8uOY8F/89oU1ZWluXtG0KuBDdCiKQS+ONPKWVVHhkxAEqpt0Ll1DPH5swCNmut\nw1rcMyECHIypaWUBwc0CjJHglk1Ns6uaRjePvtbEsUYPuVmKq8/N4MLT0hmaEd204B5a0GiUrAQu\nhBC2sra8k7XlnWDRWlo2UAaglDobY9r5Pu8OpdS/AlprvVApNUkp9a9a66f7OqBtApxQU9nMqWkr\nzW3eTMZ1Vq5/YWf7a9wca/Rw47xMLjo9nRSnNQFHFlfQyF5z5pR/xmMH4KGTT6nmAnL5HllchYMh\nljyvEEIkm1hPE79kUQaXLMrg+aeOhZ34zubylVJ7MKeeK6X2+vXoLMZMQmiO0akP54C2CXD6msoW\nLMHgYOEdU+Vya8JcGzUseXyfHG7mBL/nBE/hoQljZVpvvhsXLiqp5b+o5R5G8QJDmGVdBYQQIknY\nNQ9OAsnVWvsy0Ab01AwPWAczrCs38VjbXURo4ggnY/OdrHm/jeV/aODlf7TS0GJN0j0nBeTxfU7h\nQ3pL9qdpQdNFJ9steU4hhBAiQOASEFWcTDwY+OUU1k99CXASwPChTu67IYfvXzWUopEpvL65nbv+\n0MBftrVZ9hyKVOTtIIQQIk7OUUpdopTKMW//6rcvL6Bsr0tL+JNvtATh9kBzm4emNqPnJi1FkZFm\n3eDfTnZDyKUYZJkGIYQIxe1LshG7WxJ5APgOxuWnKow1KZVS6gFguVLqAaXURDPwCbakRA+2GYMj\nelfb5Oahl5qoa/ZQmOPgxnmZnD89nSEDDHA0blpZSwOP08VOjHjXjXeQsSEVhYN05pLJgoE1RAgh\nhAjCXGDz+iC7XgRQSm0B7gIqBroWVdJpbGyktLSURYsWsWjRonhXJyKpTsWwLEVdMzS3a2qbPDS1\nexiS5qS8vJzy8nK6uroibt8xltDGe2haA/Z4MFL+KXK4mRyWkMp4q5slhBBx4/3sxKJp1rJUgzXM\nHhrvr+lNWutHwDcm5zsRHUtbOS3HxubMmaMTOdEfwN7PXazf3s7mvZ1oYNGcDBbNNRLxZWVl0dLS\nEtHxqimmi0963Z/OFxkT3ppmQgiRkJRSm7XWc/ouGdoZc9L1morYL2czQ+23pP7xppQahrH+1BaM\nlcrBSAkzG7hUa30i0mMOmh6cZDB5VAqTR2VT1+zh56+eYOehLl+AEw2S3k8IIcJjTBOXr9QBWKK1\nXhi40VxkuxR4JNIDytlIQPnZDnIyBz4+XJFO9/E2/pzmfiGEEHbzTnkL75S3QPJkMg66aLeZ6Ldf\nC3rLLKpBrJDHyOIqFOkoM0uxIgPFELIpYTgr4lxDIYRIHLGcPXXBohz+T9logGTJZBxqvEy/xtJI\nD84glsYURvAr3NRygv+hhVfJZjFD+SrOpPlRIIQQIgFMVkpNDMhYjFJqIsZYnIhJgCNwMpw8biOP\n2+JdFSGESEiyVMPAaK0fVkqtUkpNArwLaxdhrDvZ6yrjoUiAk0C01hyqc7N9XxdH6t2MzpP/TEII\nIZKD1vp6c+HtS81ND2qtt/b3eBLg2FynS7OzuouP9nexfX8Xdc3GgOAJhU4uOTPM1b11BzT9Eeof\nAFyQdxdkfw0c0ZuBJYQQg4kGyYNjgWALayul/tObDycSEuDY2O7DXTz2WhOdLuP+GeNTuXJOBmdO\nSCU3K4zx4doNtf8HGn8JeEA3G9trbjduOaVQ8ACotKi1QQghhPCnlLoTKA63OEYuHAlwepOImYzz\nsx3MOCWVT6u76HQZAY/TARrNWRPSyM1yhM5k3LkdGh8FHbAopzfQOfFryL4KMi6KbcOEECLOrM5k\njOTBicRw4CFOjrXpy/L+PIlkMk4AXS7NrsNdbN9nXKaqbTIuUy38whAWn58J9JLJuH0LHLoYeksA\nqYbB6Jchc34Uay+EEPZlVSbjGXMy9W8rpllRpYicp7YlXCZjpdTMSMbWKKUmmUs1RETCzQSQmqI4\nY3waZ4xP46tac7jeza/fbGZfjSvk4zQujLH9QgghoinWs6g2lDewsbwREjDRX6QDh/sT3IAk+ks4\ntU0etu/rorm998BFe47jbr8fd/tlQBNaBTnNKhvwQOrE6FVWCCFEVMxblMvysgmQPIn+LCc9OAmg\nocVDxZ5ONu3poPKoG4CikU4Wnt1zFpW77TZ0VxmgwNGGuwBUuwdHqwK3BpWBco6G/B/KTCohhLCQ\n5MGxFwlwbG5bVSe/+nMzGshMV1xzbgZzp6ZRmBP8P5Lu+g3QfnKDAp0B7iEaXKBSrsaZ8xwouXAl\nhBAiecklKpsbO9zJzKJUUpzQ2qF575MO3t3RwYHjLiIaIK6AVCB9rAQ3QggxyCilipVSPRYYVEqt\nMFfsjuRYRebjSs1buFO+Y0p6cGyuMMfJrZcNpa1T88HuDl7+RxtvbW3nra3tFJ+Vzg3zssI6jga6\nnCmgjpKBRsnQYyGEsIxdl2owg48FwCyCT8suAZapnj98l2qty4IcrwhYqbVeYN5fhjGNe52V9baC\nBDg25tGaAzVudhzoYsfBLiqPunB7ID0Vpo9N5exJPRP0qdRSdNevMLps2tBAW+oQmtKy8TgcwLs4\nmEsO3yOLEhxkxrhVQgghYkVrvQ5YZ/beBOupWQOsDNgWNLgxrQgoX2bebEcCHBv7xRvNfLS/C4AU\nJ8w/PZ2ZRWlMHpVCijN4D4wz46foIf+Fp/NpXF2PcCxLA0608l7OasfNQRq4jwZ+xGjeJoUJsWmQ\nEEIkKY1KuKUazEtTK7XWlX7bSoEHQpQvAZZ4t2mtG6Jdz/4aNGNwvJmMzcyVCeGCGemcMT6VtBRw\nueHdHR28vrmN9dvbqTbH4JSXl1NaWurLZFxeXo5S+TjTl0HWW0COX3BzkqYVcOCiOubtEkKIePN+\ndpKAeWSsorVuCAhuZgGVIYIWb0LBIqVUiXnrMa7HLiSTcQLocmk+O+LyXao6VGdMFZ83I51vXGyM\nwQmWybiT7RzlWjRNQY+ryKGQZxnC+dFtgBBC2JRVmYxPnTNUP14R+4TCl6t3wqq/9xKV1nppiDIr\n+9hfAqwGFpiXvrw9PrNDPS5e5BJVAkhNUZw2LpXTxqWyGKhr9vD4a00ca3T38UgFdIXY32WWsYd2\nTrCb9ThIYSoXk052vKskhBC28kbZYd4sO+y/qcCK45rBy+Y+ijWAb1yPVwXGmBwJcMTA5Wc7yBrS\nd2CSymkM5d9o4mnAjcbo4VFkAYpsbiadudGtbBhqqeQjXqHa939LsY0/MpHzOJ1F5DE+rvUTQoi+\nxGoW1b+UjuNfSsf57l+p1h+36NBLMQYQh1IXZFsDGONz7DYeRwKcJKZwkstyhnEHrZTTyGOAy5xB\ndTWKnpmQY62OfbzBvWhcaLpfLq1kA1X8nav5OUMZEacaCiHEoFBMH70wWustSqmGgGAm19xnq+AG\nJMAZFBRpZHEdWVwX76r00EEzTlLpCnIpTeMhhXS6aI1DzYQQInx2zYMTDjO3TW/7ZgFz/KaNlwHX\nc3JqeDHwUHRr2D8S4AghhBBJygxQijGmd3sT863TWm8JUjzYJahi4AbMgEZrvdzMYrzM3D9ca73c\n+poPnAQ4Iq5SSMeDq9f9blw46ZnQUAgh7MaOeXDMQGYLIXpZzKniQQd2aq0fCnysXQOaQIMmD46w\npwKmcCG3UcAUnKShcKBw4iSVkczgUpYxjDHxrqYQQogEIz04Iq4UivHMYTxzaKCaHbyGk1RO48vk\nMCre1RNCiLAYY3DkK9VOBk0PTiJmMg5HsEzGiSqXUzif7/BFvi3BjRAiqiSTcfKTTMYJ6uGXTwBw\n59U5QPBMxgAeXDjkV4UQQgRlVSbjojl5+icV8y2oUXg2lx9hS/nnvP3U/j1a66kxe+IEMmh6cAYT\nDy6OspYP+TrvcQmf8H9pZk+8qyWEEMIisxeNZknZTIDGeNfFrmz3014pVYyxzsXygO25QKl5dy7G\ngmAJMZI7lg6yin08g8aN28wfc4z11PAumUxkBveQzeQ411IIIZJPoubBSVa2CXC8gQ0wC6gMUuRu\n/4BGKbVZKbVCgpzuqliJm7Zu2zRuNG6a2cXnvMkU/j1OtRNCCCFiwzaXqLTW68xgpUfyIbP3pjZg\n80pO9uiIsA2OMVdCCCEGN9sEOH0oAlaYGRm96jDXwBDgcvWeLC/RlZWV9V0ogUn7Epu0L+FZshq3\nRuHCGfOb6F1CBDhmJsYFAamlFwDrennIoOMNcPI5FwfpEPDGd5CBgyHkMjMOtRuYZP+AlfYlNmlf\nwiu04iDePDixvoneJUSAA8YlLO/f5iWr6wEZfxPgDH7CXJ5lNF/GQToOhpDOSKbwH8zjdQo4P6Jc\nOXYoGyk71Dla7Yv0uHZoX0ND+IsMR/M1TrTzZ5fXIlrHtUNZkdwSJsAJ8BSwuJfFwga9TCYwnbs4\nn9eYyS84jz8xlmtwkgHY44NFApz+scuXeiRlGxvDn8Vqly91O5w/u7wW0TquHcqK5Ga7RH9KqRVA\nrtZ6aYj9a/17dMI8bhPdA7oa4Hi/Kxp/BXTvWs2cPXt2yAccaz5GbkYubc1tDBsWXvLOxsbGuJet\nqamhsDD8XmQ71Dla7YvkuJGWj1bZAwcOMH78+LjWIZrHjtb5s8trkWjti/S12Lx5M1rroAtNRkIp\n9SYWjeeJ0HGt9WVxeF7701rb6gasAFb2sq8EKPa7Pyve9bXLLTMzU4dSVVelix4u0o///XG9ZMmS\nkGX92aHs7Nmzwy4bzXrYoX2RHDfS8tEqW1BQEPc6RPPY0Tp/dnktEq19kb4WQIu2wWe43Ky/Jcwl\nKjNPTj5QoZTKNcfhFMe5Wgmjss5ILXSk6QiLFi0K+3F2KBspO9Q5Wu2L9Lh2aF8kv6aj+Ron2vmz\ny2sRrePaoaxIbra5RGVOAS8GvJemVgLrtNZbzGCmPsjD1mitF8eqjnaWlZWlg61F5VX2YRkr3lvB\nhRMv5Lclv41hzQZuzpw5JNM6YoGkfYlN2pfYlFKtWuuseNdDWM82c8y0MWB4C/BQkH0NwICvkSaz\ngoLQl3731u0FjB6cWHv+n89TUV3Bz778s3493lzxN2lJ+xKbtC/h1cS7AiI6bNODIwamr9XSS54r\nYeuRrQxNH8q2726LYc3g2y9+m7/t+xsf3fYR6SnpMX1uIYQIxarVxIX9JMwYHNF/Wmv21O3BoRw0\ndTTR3Nkc0+evbqzGrd3sq98X0+cVQggxeEmAMwgcbz1OU0cTXxj9BQCOnIjdZSqtNdUnqgH4rPaz\nmD2vEEKIwU0CnEFgT+0eAC6YcAHQ9ziclz95mY+PfmzJcx9vPU67qx2Az45LgCOEECI2JMAZBPbW\nGgOM502cB4QOcFweF3e/dTcr/7HSkuc+2HjQ97c30BJCCCGiTQKcQWBP3R6y07I5c9SZKBSfN33e\na9mquio63Z3sOr7LkueubjQuT03KmySXqIQQQsSMBDiDQGVdJUX5RaQ50yjMKgzZg7Pz+E4Aquqr\naO9qH/BzewOc+UXz2Ve/j05354CPKYQQQvRFApxBYG/tXibnTwZg1NBRIQOcXTVGz41HeyzpcTnY\neJDhmcM5a9RZMpNKCCFEzEiAk+SaOpr4vPlzpgyfAsDooaND9+DU7GRo+lDgZLAzENWN1YwbNs73\n/DLQWAghRCxIgJPkvGtQFeUXAScDnN4SPO6s2clFky4iIyXDd7lqIKobqxmbM5aivCIcyiEDjYUQ\nQsSEBDhJzrtEg38PTmtXK00dTT3KNrY3cqTpCKcVnsbUgqkD7sFxe9wcbjrMuGHjGJI6hHHDxslA\nYyGEEDEhAU6S21u7l1RHKuNzxwNGgAPBp4rvPr4bgGmF05heOJ2dNTu79fS4PC6e2vQU9W3B1j3t\n6WjzUVweF6cMOwWAqQVTpQdHCCFETEiAk+T21O1hYt5EUhzGuqqjho4C4HDT4R5ld9YYl6SmF05n\nWuE06trqON563Ld/4/6NPPjug6z+aHVYz+3NgTNu2DgApg6fSlV9FV3urv43SADQ2tmKR3viXQ0h\nhLAtCXCS3N7avb7xNxC6B2dnzU5yh+QyMnsk0wum+7Z5rd2zFoCKQ70v6unPG+B4e3CmDJ+Cy+OK\naCZVp7uTQ42H2Hp4a8j8PYNJS2cLF5RdwB+2/iHeVRFCCNuSACeJdbo7OdBwgMnDJ/u2jcgegUM5\nggYLu2p2Ma1wGkopphVOA04GOB7tYd2edQBsqt4UVu9BdWM1CsWYnDGA0YMD4WU0bupo4qrfX8WM\nn8/gwqcupOT5Em5Zc0ufjxsMPjz4IQ3tDbz12VvxrooQQtiWBDhJbH/9ftzazZT8Kb5tKY4URmSP\n6NGD49Eedh3fxbQCI7DJy8hjVPYo30Dj7Z9vp6alhgsmXsCJjhO+8TqhHGw8yKiho0hzpgHGTC6F\nCmug8Ys7XmTHsR2Uzi3lJwt/wk1n38RntZ+xv2F/2O1PVhsPbARg6+GttHa2xrk2QghhTxLgJDHv\nDCr/S1QQPBfOwcaDtHa1Mr1wum/bqQWn+pZsWPvZWlIcKSy7cBlg9CL0xZsDxysjNYNxuX3PpPJo\nD/+z9X+YOXomyy9azg1n3cAts28B4N3Kd/t83mS3Yd8GctJz6HR3sunQpnhXRwghbEkCnCT21FkJ\nswAADuZJREFUYfWHpDnTfFmMvUZn9wxw/AcYe00vnM6e2j10ubtYu2ct55xyDqeNOI0xQ8eE9cVa\nfcLIgeNv6vCpfSb727BvA1X1Vdw08ybftol5E5mYN5F3qt7p83mT2dHmo3xW+xnfmv0t0pxpbNy3\nMd5VEkIIW5IAJ0l1uDp49dNXKZ5STGZaZrd9wZL97a7ZjUL5xsmAEeB0ujt5e+/b7K3by4IpCwCY\ne8pcNlVv6jVZoPf5jzYd7daDA8ZA4331+0LOpPrD1j9QkFnA5dMu77Z9/qT5fHDwA0vWyIqlAw0H\n+OmGn9LW1TbgY23cbwQ0l065lNljZ7Nh/4YBH1MIIZKRBDhJ6u29b1PfVs/iMxb32Dc6ZzQdrg4a\n2ht823bW7GRC3oRuwZC3N+cXH/wCgOIpxYAR4NS01LCvYV+vz3+46TAa7ZtB5TV1+FS6PF38acef\nOHzicI8g6UDDAf5a+VduPOtG39gdr/lF8+lwdfDBwQ/6fgFsoqalhm+s/ga/+uBXPP/P5wd8vI37\nN5Kfkc/0wunMmzCPXcd3UdNSY0FNhRAiuUiAk6RWf7yaUUNHcf6E83vsCzZVfOfxnb4Bxl6T8ieR\n6kjlk2OfcMbIM3yzoc455RzAmE3VG+8q4oE9OGePOZshKUP44V9+yAVlF3Duk+dy37r7ON5i5Nt5\nbttzOB1Ovnb213oc85xTziEjJcOSy1RHm4/yuy2/w+VxDfhYvWnqaOLbL36bmpYaTi04lac3PU2H\nq6Pfx9Nas3H/Rr404Us4lMN3br29OkIIIU6SACcJfd70OX/b9zeuO/06nA5nj/2BAU5rZyv76/d3\nG38DGON3zCnm3t4bMAYt52fkhxxoHJgDx2tS3iQ2//tm1nxtDfddch/nTzifF7a/wCVPX8IT7z/B\n6o9Xs3DqQkZmj+xxzPSUdL404Uu8W/luyMtjffFoD7e9dhs/fvvHrPxwZciyLo+L96re6zM5YZe7\ni2cqnuHZLc9SUV1BfVs9t75yK7uO7+KXX/kl91x8D8dajvHijhf7Xe/dx3dT01LDvAnzADh95Onk\nZeSxYZ9cphJCiEAp8a6AsN5Ln7yER3u47ozrgu73BTgnjABn2+fb0Ghf7ht/0wqmsbNmJwunLPRt\nU0pxzinnhEz4V91YTaojNWigMiR1CDPHzGTmmJnczM1877zvseK9FTy68VEAbp55c6/HnT9pPuv3\nrqeqvqrH7LBwPbv5WTZVb6Iov4gn/v4ElxRdwowRM3qU01pzz1/uYfXHq7nmtGt4+PKHUUoFLXfv\nuntZ9dGqHvt+esVPuWjSRWit+cKoL1D2YRnXn3m9L7O0vy2Ht3DPX+7hP877D66YdkWP/d7xNt6e\nG4dycN7489i4fyNa66B1E0KIwUp6cJKM1prVH63m3FPOZULuhKBlhmcOJ8WRwpGmI2w7so3vvvpd\nCrMKfZee/F17+rXceNaNnFpwarftc8fN5WDjQQ6f6LnkAxgBzpicMUF7kAJNyp/Er6/+NS/c+AI/\nuvRHzBk7p9eyF026CIB3Kt/p87jBVNZV8siGR7h08qWs+uoqhg0Zxp1/vpNOd2ePso///XFWf7ya\nmaNn8tInL/GzjT8LeszfbP4Nqz5axa3n3sr733mfp695mtvPv50nv/IkV592NWAEhf/2xX/jYONB\nXt/5eo9jVByq4JbVt/BZ7Wfc/trtvLbztR5lNuzbQFF+ke9SIcC8CfM41nJMFjEVQogAEuAkiZoa\nY6BpxaEK9jfsp+TMkl7LOh1ORmaP5K+Vf+WmVTeRk57Dqq+uIi8jr0fZeRPncf/C+3v0DswdO9f3\nfMEE5sAJx9xT5vL1mV/v8VxlZWW+v8cOG8vU4VP7NQ7H5XFx55/vJCMlg/sX3k9eRh73L7yfT2s+\n5Zfv/7Jb2T9u/yOPv/84JWeUsPprq7nhrBuCDhRev3c9D7zzAJdNvYzvz/s+I7JHcPHki/nued9l\n4dSF3cpeMvkSTi04lSf/8WS3TND3/upevrXmWxRmF/LWN99i1thZ3PH6Hbz66au+Mh2uDj6s/rDH\nmCrv5So7j8PxP3/JSNqX8AriXQERHXKJKkkcaz3Gj9/+MVsPbyU7LZvLpl4WsvzooaOpOFTBtIJp\nPFvyLCOyR0T0fNMLp5Odls0zFc+w7ci2Hvv31O7hyhlXRnTM3pSVlVFaWuq7f9Gki/jdlt/xo/U/\niuiyzKETh9h2ZBuPXvkohVmFgDG26JrTruHJfzxJQ3sDToeTTncnq7av4sKJF/LfC/4bpRQ/Lv4x\nR5uOct+6+9hxdAfpKelorVnz8RpOH3k6j1zxCA4V+veCQzm49dxbueP1O/jBGz8gLyMPj/bwXONz\nTBoxiedueI6R2SN55tpnWPKnJfzgjR/w/oH3yUjNoL6tnnZXuy+g8Ro7bCwTcifwv//8X9+4J7t5\nYdMLfD4ledcRk/YlvMJ4V0BEhwQ4SUKnaV7a8RIAt8y6pUfuG3/l5eVcOOlCMlIzeOzKxxg2ZFjI\nsosWLeqx3elwcs3p1/DKJ69woOEAAF1dXaSmpgKQ6kzt8WUcznHDcdWMq3jl01d4+ZOXg+73r0eg\nm86+iSunnQy8ysvLuXfBvVTVV3XrMfni+C/yxFVPkOo0jpPiSOHylMtxT3Dzxq43fOUm5U2i7Joy\nMlIzwmrfFdOu4Lltz3W7xOZsdPL8d5/3BZlZaVk8fe3T3PH6Hby5+01fucn5k6nfUQ9Tuh/zhjNv\n4Ncf/tp3/sN9LWJVtmNiR9C6xbIO0Tx256TOuLcvmq9ForUv0tdCJDGttdyS4JaZmanDtWTJkoQq\nO3v27LDLRrMedmhfJMeNtHy0yhYUFMS9DtE8drTOn11ei0RrX6SvBdCibfAZLjfrb8o4vyLRKaWa\ngF1hFh8GNCZQ2QLgeJhlo1kPO7QvkuNGWj5aZccDB+Jch2geO1rnzy6vRaK1L9LXYprWemgE5UWC\nkABHCCGEEElHZlEJIYQQIulIgCOEEEKIpCMBjhBCCCGSjgQ4QgghhEg6EuAIIYQQIulIoj8Rc0qp\nZUAlkA+gtQ6ZCz6c8kqpYmCB1nq55RUWg0o03p9+ZVdrrRdbV1shRG9kmrhNWf0hG+nxokUptRJY\nrbVe518vrfWa/pT3BjbALHP70hg0IyQrz51SKhfwrlMxF6ONcQ3iotS+BmCyuT9u7bP6/RlQdgVQ\nrLWeHbUGhMHi81cCFAFlWusGpVQpUKG13hKl6vcpCp+dRcBSYK+5qdJ7voXNxTvToNx63oCVGB+E\n3vvLgJL+lo/0eFFumw64PwvYPNDywApgZRKeuxUB5TcHbkvC9i2LY/t0wH2r3p9F5nu012Ml6Pkr\nBbTfLS6fK1FsXxGwNmD/WqvqK7covx/iXQG5BTkpFn/IRnq8KLZrFlAfsC03sH79KW+jAMeyc2e2\ndVnA/tLA1yRR22fer/dvI7Aao0ckHm2L5vuzFCi2QYBj9fkrNYOAoni2K4rtWx0Q8OQCufFup9zC\nu8kgY5tRSs3C6K73V4nxHy/i8pEeL8qKgLpgO8xLFQMtH1dWnzvMX/1mOa86jA/ZmItC+wBmA/6X\nBIqATQOrab9F5f1pXkZdZUUFByJK5w+tdaXWutKqevZXFD47c4ESwHc5SmvdoLUOfIywKRlkbD8h\nPzSD/Ofq60M20uPFSz49P2ysLB8Llp47rfUWpdQC3X08wwL8PnBjzOr2Nfh/MZpfOLla64esqrCF\nBvL+zNXG+BTraxUZy8+f+XepWS6f+J4/qz8753jLmeNwAOZqmciQMKQHJ7HkW1w+0uNFU9APGgvL\nx1u/zp32G8xofuheD9jxA7bf702lVK75JbkCsOsMo369P5VSJbqXAco209/zt05rXaa1XqONwbnD\nzUG7dtOf9nl74fLN9q0B9poDy0UCkAAnsVgdBMQ6SPDNVPDj/SIP9us40vJ2ZsW5ewpYrOM4QyWE\nfrfP7Mkp01ovAJ4yg514sPT9af7qT5T3ab/OX5BLU5uAuy2pkbX6074G6P4jA6jg5KxGYXMS4NiP\n1UGAbYIE84s5cCxDLhD0CzvS8jYQtQDOnGK8Usd3eqrl7fPr+vdaad5iLgrvz1nALKXUMrNXYymQ\na94PbHcsWB7AKaV0QFsaiNMYMax/fwYLirzvW9uNARQ9SYBjM1Z/yNowSCgzc2d43YDfF5pSalbA\nL/iQ5e0kWgGc2f61+mSulXgMELe8febg273Bvuzj+AVi2fvTvKzxkPcGrAUazPsxH5QbhfdnHbA8\noC25GIFDzEXps7Mh4L2Ya+5LlJ65QU0CHHuyOgiwTZCgjUR8c5VSJWYbanX3xFrFGL90wypvvhbL\nMGY7FJu/juMSAJgsPXdmEJAPVJhjVXIxXqN4sbJ9FRgJ4vy/EBdgjOuIyxeI1e9PL3PfYowBq8uS\nJIALdo6WEt8xYpZ/dmKMe/MqBuw4CF4EIZmMbcq8JLGJIDMTzC/0G7RfRtRQ5cPZL6xj1bkzvwTr\ngzzFGh3HdP9WvjfNYLSYk5mMc7UNslEnM4vPX2Am6r3BArxYitJnZ615d7jMokocEuAIIYQQIunI\nJSohhBBCJB0JcIQQQgiRdCTAEUIIIUTSkQBHCCGEEElHAhwhhBBCJB0JcIQQQgiRdCTAEUIIIUTS\nkQBHCCGEEElHAhwhhBBCJB0JcIQQQgiRdCTAEUIIIUTSkQBHCCGEEElHAhwhhBBCJB0JcIQQQgiR\ndCTAEUIIIUTSkQBHCCGEEElHAhwhRNQopdYqpbRSqsS8X6KUyo13vYQQyU9preNdByFEElJKLQO2\nAHVAMTAc2KS1XhPXigkhBgUJcIQQUaeUmgUUSXAjhIgVCXCEEFFlBjf5Wut18a6LEGLwkDE4Qoio\nUUoVA0hwI4SINQlwhBBRYQ4srtNabwnYJoQQUZcS7woIIZKP2XMzFyhSSnk33w0sj1ulhBCDigQ4\nQghLmdPAi7TWy5VSRcBmjJlUS7XWlfGtnRBisJBBxkIIIYRIOjIGRwghhBBJRwIcIYQQQiQdCXCE\nEEIIkXQkwBFCCCFE0pEARwghhBBJRwIcIYQQQiQdCXCEEEIIkXQkwBFCCCFE0pEARwghhBBJ5/8D\nS/TC+ndevy4AAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fafcc74bcd0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Create a GridSpec\n",
"\n",
"gs = grd.GridSpec(2, 4, height_ratios=[10,1.5], width_ratios=[.75,6,.25, .4], wspace=0., hspace=0.)\n",
"\n",
"# Image plot\n",
"\n",
"ax1 = plt.subplot(gs[0,1])\n",
"\n",
"ZD = np.log10(J_heatmap/(GeV**2*Centimeter**-5))\n",
"x = X.ravel()\n",
"y = Y.ravel()\n",
"z = ZD.ravel()\n",
"\n",
"gridsize=30\n",
"\n",
"p = plt.hexbin(x, y, C=z, gridsize=gridsize, bins=None, cmap = CM.jet)\n",
"plt.axis([x.min(), x.max(), y.min(), y.max()])\n",
"plt.axis([0,0.06, 12,16.5])\n",
"ax1.xaxis.set_ticklabels([])\n",
"ax1.yaxis.set_ticklabels([])\n",
"\n",
"ax1.text(0.002, 15.85, r\"$\\mathbf{J > 10^{17.5}$ \\textbf{GeV}$^2$ \\textbf{s}$^{-5}$\", fontsize=18)\n",
"\n",
"# Color bar\n",
"\n",
"colorAx = plt.subplot(gs[0,3])\n",
"cb = plt.colorbar(p, cax = colorAx)\n",
"cb.set_label('log$_{10}(J_\\mathrm{mean}~[\\mathrm{GeV}^{\\mathbf 2} \\mathrm{s}^{\\mathbf {-5}}])$')\n",
"\n",
"cb.set_clim(16.,18.8)\n",
"\n",
"# z profile plot\n",
"\n",
"ax2 = plt.subplot(gs[1,1])\n",
"ax2.yaxis.set_ticklabels([])\n",
"ax2.set_yticks([])\n",
"ax2.plot(z_bins, np.sum(J_heatmap, axis=0)/(GeV**2*Centimeter**-5),color='forestgreen')\n",
"plt.xlim(0,0.06)\n",
"plt.xlabel(r\"$z$\")\n",
"\n",
"# Mvir profile plot\n",
"\n",
"ax3 = plt.subplot(gs[0,0])\n",
"ax3.xaxis.set_ticklabels([])\n",
"ax3.plot(np.sum(J_heatmap, axis=1)/(GeV**2*Centimeter**-5), np.log10(M_bins), color='cornflowerblue')\n",
"plt.ylim(12,16.5)\n",
"ax3.set_xticks([])\n",
"plt.ylabel(r\"log$_{10}(M_\\mathrm{vir}~[M_\\odot])$\")\n",
"plt.gca().invert_xaxis()\n",
"\n",
"plt.tight_layout()\n",
"plt.savefig(\"DarkSky_heatmap_Jcut.pdf\")"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.7"
},
"latex_envs": {
"bibliofile": "biblio.bib",
"cite_by": "apalike",
"current_citInitial": 1,
"eqLabelWithNumbers": true,
"eqNumInitial": 0
}
},
"nbformat": 4,
"nbformat_minor": 1
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment