{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Going beyond galaxies as tracers with `halomod`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`halomod` is written in a way that is most native to applications of halo models of galaxies. Therefore, modifications and extensions in the context of galaxy clustering (as well as HI assuming HI is trivially related to galaxies) are very straightforward. However, it may not be as straightforward when dealing with other tracers. In this tutorial, we use the flux density power spectrum of [arxiv:0906.3020](https://arxiv.org/abs/0906.3020) to demonstrate how to fully utilise the flexibility of `halomod`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The flux density power spectrum can modelled as (see Sec 2.3 of [arxiv:0906.3020](https://arxiv.org/abs/0906.3020)):\n", "$$\n", "P_{1h}(k) = |u_J(k)|^2 \\int_{M_{\\rm min}}^{\\infty} {\\rm d}m\\, n(m) \\bigg(\\frac{m}{\\bar{\\rho}_{\\rm gal}}\\bigg)^2\n", "$$\n", "\n", "$$\n", "P_{2h}(k)=|u_J(k)|^2\\bigg[\\int_{M_{\\rm min}}^{\\infty}{\\rm d}m\\,n(m)b(m)\\Big(\\frac{m}{\\bar{\\rho}_{\\rm gal}}\\Big)\\bigg]^2 P_{\\rm lin}(k)\n", "$$\n", "\n", "where $u_J(k)={\\rm arctan}(k\\lambda_{\\rm mfp})/(k\\lambda_{\\rm mfp})$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## HOD\n", "Once we have the expression of the power spectrum we want, we should try to identify the halo model components. Comparing it to the standard halo model formalism, it's easy to see that it effectively means:\n", "$$\n", "\\langle M_{\\rm cen}\\rangle \\equiv 0\n", "$$\n", "\n", "$$\n", "\\langle M_{\\rm sat}\\rangle \\equiv A_{\\rm sat}\n", "$$\n", "\n", "where $A_{\\rm sat}$ is a constant so that the total satellite occupation is equal to the mean mass density of galaxies:\n", "\n", "$$\n", "\\int_{M_{\\rm min}} {\\rm d}m\\, n(m)A = \\bar{\\rho}_{\\rm gal}\n", "$$\n", "\n", "This HOD has already been defined within `halomod` by the `Constant` HOD class:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "import hmf\n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", "\n", "import halomod\n", "from halomod import TracerHaloModel" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using halomod v2.0.2.dev51+geee0902 and hmf v3.5.0\n" ] } ], "source": [ "print(f\"Using halomod v{halomod.__version__} and hmf v{hmf.__version__}\")" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "hm = TracerHaloModel(hod_model=\"Constant\", transfer_model=\"EH\")" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hm.central_occupation" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAh+UlEQVR4nO3df3BU1f3/8dcmkA0oCdrIhsBK0IqoQMAgaaB+/bUlIoMw/WFExyBFrHzRAVOnEJVEamv8BWWqkSgFseMoqKPYET7wgRS0SpSakE/VKooiQWEX0C/ZGCSh2fv9w2H9xPwgd3PD2d08HzM7Y27Oufs+Hnfz8u65Z12WZVkCAAAwJMF0AQAAoGcjjAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwqpfpAjojFApp//796tevn1wul+lyAABAJ1iWpfr6emVkZCghof3rHzERRvbv3y+v12u6DAAAEIF9+/Zp8ODB7f4+JsJIv379JH03mJSUFMPVAACAzggGg/J6veG/4+2JiTBy4qOZlJQUwggAADHmZEssWMAKAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjLIdRt544w1NmTJFGRkZcrlcWrdu3Un7bNu2TRdffLHcbrd+/OMfa/Xq1RGUCgAA4pHtMNLQ0KCsrCyVlZV1qv2ePXs0efJkXXHFFaqpqdH8+fN1yy23aNOmTbaLBQAA8cf2d9NMmjRJkyZN6nT78vJyDR06VEuWLJEkXXDBBXrzzTf1pz/9SXl5eXafHgAAxJlu/6K8yspK+Xy+Fsfy8vI0f/78dvs0NjaqsbEx/HMwGOyu8gDAEe99Uad1NV8qZFmmSwEi8usJQ+U9s6+R5+72MOL3++XxeFoc83g8CgaD+vbbb9WnT59WfUpLS7V48eLuLg0AHFPyt/dVXXvEdBlAxKZkZcRvGIlEUVGRCgsLwz8Hg0F5vV6DFQFAx442NUuSrs3KkPfM1v+TBUQ7T0qysefu9jCSnp6uQCDQ4lggEFBKSkqbV0Ukye12y+12d3dpAOC468Z69dPz0kyXAcSUbt9nJDc3VxUVFS2Obd68Wbm5ud391AAAIAbYDiPffPONampqVFNTI+m7W3drampUW1sr6buPWAoKCsLtb7vtNn322Wf63e9+p48++khPPPGEXnjhBd15553OjAAAAMQ022Hk3Xff1ZgxYzRmzBhJUmFhocaMGaPi4mJJ0oEDB8LBRJKGDh2q9evXa/PmzcrKytKSJUv0l7/8hdt6AcQVbqIBImd7zcjll18uq4NXXVu7q15++eXauXOn3acCAAA9AN9NAwAOcrlMVwDEHsIIAAAwijACAACMIowAgAMssYIViBRhBAAAGEUYAQAHsX4VsI8wAgAAjCKMAAAAowgjAADAKMIIADiA7eCByBFGAACAUYQRAHASt9MAthFGAACAUYQRAABgFGEEABzA+lUgcoQRAABgFGEEABzkYgUrYBthBAAAGEUYAQAARhFGAMABFluwAhEjjAAAAKMIIwDgIBfrVwHbCCMAAMAowggAADCKMAIADmD5KhA5wggAADCKMAIADmL9KmAfYQQAABhFGAEAAEYRRgAAgFGEEQBwArfTABEjjAAAAKMIIwDgIBf7wQO2EUYAAIBRhBEAAGAUYQQAHMD6VSByhBEAAGAUYQQAHMT6VcA+wggAADCKMAIAAIwijACAAyyLJaxApAgjAADAKMIIADiI9auAfYQRAABgFGEEAAAYRRgBAAewfBWIHGEEAAAYRRgBAAexAytgH2EEAAAYRRgBAABGEUYAAIBRhBEAcAC7wQORI4wAAACjCCMA4ChupwHsIowAAACjCCMAAMAowggAOMBiQ3ggYhGFkbKyMmVmZio5OVk5OTnasWNHh+2XLVum888/X3369JHX69Wdd96pY8eORVQwAACIL7bDyNq1a1VYWKiSkhJVV1crKytLeXl5OnjwYJvtn3vuOS1cuFAlJSX68MMPtXLlSq1du1Z33313l4sHgGjDdvCAfbbDyNKlSzV79mzNnDlTF154ocrLy9W3b1+tWrWqzfbbt2/XhAkTdMMNNygzM1MTJ07U9OnTT3o1BQAA9Ay2wkhTU5Oqqqrk8/m+P0FCgnw+nyorK9vsM378eFVVVYXDx2effaYNGzbommuuafd5GhsbFQwGWzwAAEB86mWn8eHDh9Xc3CyPx9PiuMfj0UcffdRmnxtuuEGHDx/WT3/6U1mWpf/85z+67bbbOvyYprS0VIsXL7ZTGgAYxQ6sQOS6/W6abdu26YEHHtATTzyh6upqvfzyy1q/fr3uv//+dvsUFRWprq4u/Ni3b193lwkAAAyxdWUkLS1NiYmJCgQCLY4HAgGlp6e32WfRokW66aabdMstt0iSRo4cqYaGBt1666265557lJDQOg+53W653W47pQFAVGD9KmCfrSsjSUlJys7OVkVFRfhYKBRSRUWFcnNz2+xz9OjRVoEjMTFRkmRxXRMAgB7P1pURSSosLNSMGTM0duxYjRs3TsuWLVNDQ4NmzpwpSSooKNCgQYNUWloqSZoyZYqWLl2qMWPGKCcnR7t379aiRYs0ZcqUcCgBAAA9l+0wkp+fr0OHDqm4uFh+v1+jR4/Wxo0bw4taa2trW1wJuffee+VyuXTvvffqyy+/1FlnnaUpU6boj3/8o3OjAADDuNALRM5lxcBnJcFgUKmpqaqrq1NKSorpcgCglQkP/l1fHvlWr/zf8Rpz9hmmywGiQmf/fvPdNADgIBdbsAK2EUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQBwEPfSAPYRRgAAgFGEEQAAYBRhBAAcEAObWQNRizACAACMIowAgIPYDR6wjzACAACMIowAAACjCCMA4ACWrwKRI4wAAACjCCMA4CAXe7ACthFGAACAUYQRAABgFGEEABzABqxA5AgjAADAKMIIADiIHVgB+wgjAADAKMIIAAAwijACAA6w2IMViBhhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAAewHTwQOcIIAAAwijACAA5iO3jAPsIIAAAwijACAACMIowAgANYvwpEjjACAACMIowAgINcYgUrYBdhBAAAGEUYAQAARhFGAMAB7MAKRI4wAgAAjCKMAICD2IEVsI8wAgAAjCKMAAAAowgjAOAIVrACkSKMAAAAowgjAOAgFrAC9hFGAACAUYQRAABgFGEEAAAYRRgBAAewHTwQOcIIAAAwijACAA5yidtpALsIIwAAwCjCCAAAMCqiMFJWVqbMzEwlJycrJydHO3bs6LD9kSNHNHfuXA0cOFBut1vDhg3Thg0bIioYAKIR61eByPWy22Ht2rUqLCxUeXm5cnJytGzZMuXl5WnXrl0aMGBAq/ZNTU362c9+pgEDBuill17SoEGDtHfvXvXv39+J+gEAQIyzHUaWLl2q2bNna+bMmZKk8vJyrV+/XqtWrdLChQtbtV+1apW+/vprbd++Xb1795YkZWZmdq1qAIhSbAcP2GfrY5qmpiZVVVXJ5/N9f4KEBPl8PlVWVrbZ529/+5tyc3M1d+5ceTwejRgxQg888ICam5vbfZ7GxkYFg8EWDwAAEJ9shZHDhw+rublZHo+nxXGPxyO/399mn88++0wvvfSSmpubtWHDBi1atEhLlizRH/7wh3afp7S0VKmpqeGH1+u1UyYAAIgh3X43TSgU0oABA/TUU08pOztb+fn5uueee1ReXt5un6KiItXV1YUf+/bt6+4yAaBLLLZgBSJma81IWlqaEhMTFQgEWhwPBAJKT09vs8/AgQPVu3dvJSYmho9dcMEF8vv9ampqUlJSUqs+brdbbrfbTmkAACBG2boykpSUpOzsbFVUVISPhUIhVVRUKDc3t80+EyZM0O7duxUKhcLHPv74Yw0cOLDNIAIAsYz1q4B9tj+mKSws1IoVK/TMM8/oww8/1Jw5c9TQ0BC+u6agoEBFRUXh9nPmzNHXX3+tefPm6eOPP9b69ev1wAMPaO7cuc6NAgAAxCzbt/bm5+fr0KFDKi4ult/v1+jRo7Vx48bwotba2lolJHyfcbxerzZt2qQ777xTo0aN0qBBgzRv3jwtWLDAuVEAAICY5bJiYNVVMBhUamqq6urqlJKSYrocAGhl9O//W0eOHtfmO/+PzvP0M10OEBU6+/eb76YBAABGEUYAwEHswArYRxgBAABGEUYAAIBRhBEAAGAUYQQAHBD99yUC0YswAgAAjCKMAICjuJ0GsIswAgAAjCKMAAAAowgjAOCAGPhmDSBqEUYAAIBRhBEAcBDbwQP2EUYAAIBRhBEAAGAUYQQAHMDyVSByhBEAAGAUYQQAHMT6VcA+wggAADCKMAIAAIwijACAE1jBCkSMMAIAAIwijACAg1xswQrYRhgBAABGEUYAAIBRhBEAcADrV4HIEUYAAIBRhBEAAGAUYQQAHMS9NIB9hBEAAGAUYQQAABhFGAEAB1gW99MAkSKMAAAAowgjAOAgdoMH7COMAAAAowgjAADAKMIIADiA5atA5AgjAADAKMIIADjIxR6sgG2EEQAAYBRhBAAAGEUYAQAHsAErEDnCCAAAMIowAgAOYgdWwD7CCAAAMIowAgAAjCKMAIADLPZgBSJGGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAMABbAcPRI4wAgAAjCKMAICD2A4esI8wAgAAjIoojJSVlSkzM1PJycnKycnRjh07OtVvzZo1crlcmjZtWiRPCwAA4pDtMLJ27VoVFhaqpKRE1dXVysrKUl5eng4ePNhhv88//1x33XWXLr300oiLBYBoxfpVIHK2w8jSpUs1e/ZszZw5UxdeeKHKy8vVt29frVq1qt0+zc3NuvHGG7V48WKdc845XSoYAADEF1thpKmpSVVVVfL5fN+fICFBPp9PlZWV7fb7/e9/rwEDBmjWrFmdep7GxkYFg8EWDwCIBS5WsAK22Qojhw8fVnNzszweT4vjHo9Hfr+/zT5vvvmmVq5cqRUrVnT6eUpLS5Wamhp+eL1eO2UCAIAY0q1309TX1+umm27SihUrlJaW1ul+RUVFqqurCz/27dvXjVUCAACTetlpnJaWpsTERAUCgRbHA4GA0tPTW7X/9NNP9fnnn2vKlCnhY6FQ6Lsn7tVLu3bt0rnnntuqn9vtltvttlMaAJjFClYgYraujCQlJSk7O1sVFRXhY6FQSBUVFcrNzW3Vfvjw4XrvvfdUU1MTflx77bW64oorVFNTw8cvAADA3pURSSosLNSMGTM0duxYjRs3TsuWLVNDQ4NmzpwpSSooKNCgQYNUWlqq5ORkjRgxokX//v37S1Kr4wAQD1i+CthnO4zk5+fr0KFDKi4ult/v1+jRo7Vx48bwotba2lolJLCxKwAA6ByXZUX/d00Gg0Glpqaqrq5OKSkppssBgFaG3fNfamoOafvCK5XRv4/pcoCo0Nm/31zCAAAHWKxgBSJGGAEAAEYRRgDAQWzACthHGAEAAEYRRgAAgFGEEQAAYBRhBAAcEP2bJADRizACAACMIowAgINcbAgP2EYYAQAARhFGAACAUYQRAHAA61eByBFGAACAUYQRAHAQ28ED9hFGAACAUYQRAABgFGEEABxgsQUrEDHCCAAAMIowAgAOYv0qYB9hBAAAGEUYAQAARhFGAMABLF8FIkcYAQAARhFGAMBJrGAFbCOMAAAAowgjAADAKMIIADiADViByBFGAACAUYQRAABgFGEEABzk4nYawDbCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIADjIxfpVwDbCCAAAMIowAgAAjCKMAEAXWewFD3QJYQQAABhFGAEAB7F+FbCPMAIAAIwijAAAAKMIIwDQRaxfBbqGMAIAAIwijACAg1xswQrYRhgBAABGEUYAAIBRhBEA6CLWrwJdQxgBAABGEUYAwEEsXwXsI4wAAACjCCMAAMAowggAADCKMAIAXWSxHzzQJYQRAABgFGEEABzEbvCAfRGFkbKyMmVmZio5OVk5OTnasWNHu21XrFihSy+9VGeccYbOOOMM+Xy+DtsDAICexXYYWbt2rQoLC1VSUqLq6mplZWUpLy9PBw8ebLP9tm3bNH36dG3dulWVlZXyer2aOHGivvzyyy4XDwAAYp/LsrnyKicnR5dccokef/xxSVIoFJLX69Udd9yhhQsXnrR/c3OzzjjjDD3++OMqKCjo1HMGg0Glpqaqrq5OKSkpdsoFgG53vDmk8+75L0lSTfHP1L9vkuGKgOjQ2b/ftq6MNDU1qaqqSj6f7/sTJCTI5/OpsrKyU+c4evSojh8/rjPPPLPdNo2NjQoGgy0eAAAgPtkKI4cPH1Zzc7M8Hk+L4x6PR36/v1PnWLBggTIyMloEmh8qLS1Vampq+OH1eu2UCQDGuNgQHrDtlN5N8+CDD2rNmjV65ZVXlJyc3G67oqIi1dXVhR/79u07hVUCAIBTqZedxmlpaUpMTFQgEGhxPBAIKD09vcO+jz76qB588EFt2bJFo0aN6rCt2+2W2+22UxoAAIhRtq6MJCUlKTs7WxUVFeFjoVBIFRUVys3Nbbffww8/rPvvv18bN27U2LFjI68WAKIQG7ACXWPryogkFRYWasaMGRo7dqzGjRunZcuWqaGhQTNnzpQkFRQUaNCgQSotLZUkPfTQQyouLtZzzz2nzMzM8NqS008/XaeffrqDQwEAALHIdhjJz8/XoUOHVFxcLL/fr9GjR2vjxo3hRa21tbVKSPj+gsvy5cvV1NSkX/7yly3OU1JSovvuu69r1QNAtGH9KmCb7X1GTGCfEQDRrOk/IQ2797t9Rv6nZKJS+/Q2XBEQHbplnxEAAACnEUYAoIssRf0FZiCqEUYAAIBRhBEAcJCLBayAbYQRAABgFGEEAAAYRRgBAABGEUYAoIuif7cmILoRRgAAgFGEEQBwEDfTAPYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAMBBLvaDB2wjjAAAAKMIIwAAwCjCCAB0ETuwAl1DGAEAAEYRRgDAQSxfBewjjAAAAKMIIwAAwCjCCAB0kSVWsAJdQRgBAABGEUYAwEFswArYRxgBAABGEUYAAIBRhBEAAGAUYQQAuojt4IGuIYwAAACjCCMA4CAXG8IDthFGAACAUYQRAABgFGEEALqI9atA1xBGAACAUYQRAHAQ28ED9hFGAACAUYQRAABgFGEEALrIYgtWoEsIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjANBFLF8FuoYwAgAAjCKMAICD2IEVsI8wAgAAjCKMAAAAowgjANBFbMAKdA1hBAAAGEUYAQAARhFGAMBBLnE7DWAXYQQAABhFGAEAAEYRRgCgq7ibBugSwggAADAqojBSVlamzMxMJScnKycnRzt27Oiw/Ysvvqjhw4crOTlZI0eO1IYNGyIqFgCiHdvBA/bZDiNr165VYWGhSkpKVF1draysLOXl5engwYNttt++fbumT5+uWbNmaefOnZo2bZqmTZum999/v8vFAwCA2OeyLHt7B+bk5OiSSy7R448/LkkKhULyer264447tHDhwlbt8/Pz1dDQoNdeey187Cc/+YlGjx6t8vLyTj1nMBhUamqq6urqlJKSYqdcAOh2dUePK+v3/y1J+uSPk9Q7kU/AAanzf7972TlpU1OTqqqqVFRUFD6WkJAgn8+nysrKNvtUVlaqsLCwxbG8vDytW7eu3edpbGxUY2Nj+OdgMGinzE5b+eYeffH/jnbLuQH0HI3/CZkuAYhptsLI4cOH1dzcLI/H0+K4x+PRRx991GYfv9/fZnu/39/u85SWlmrx4sV2SovI+n/tV3XtkW5/HgA9Q9+kRCWwaASwzVYYOVWKiopaXE0JBoPyer2OP88vsgcr99wfOX5eAD3TT875kRITCCOAXbbCSFpamhITExUIBFocDwQCSk9Pb7NPenq6rfaS5Ha75Xa77ZQWkRtzhnT7cwAAgI7ZWmWVlJSk7OxsVVRUhI+FQiFVVFQoNze3zT65ubkt2kvS5s2b220PAAB6Ftsf0xQWFmrGjBkaO3asxo0bp2XLlqmhoUEzZ86UJBUUFGjQoEEqLS2VJM2bN0+XXXaZlixZosmTJ2vNmjV699139dRTTzk7EgAAEJNsh5H8/HwdOnRIxcXF8vv9Gj16tDZu3BhepFpbW6uEhO8vuIwfP17PPfec7r33Xt19990677zztG7dOo0YMcK5UQAAgJhle58RE9hnBACA2NPZv9/szAMAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMsr0dvAknNokNBoOGKwEAAJ114u/2yTZ7j4kwUl9fL0nyer2GKwEAAHbV19crNTW13d/HxHfThEIh7d+/X/369ZPL5XLsvMFgUF6vV/v27etx33nD2Hve2HvquCXG3hPH3lPHLUXX2C3LUn19vTIyMlp8ie4PxcSVkYSEBA0ePLjbzp+SkmJ8wkxh7D1v7D113BJj74lj76njlqJn7B1dETmBBawAAMAowggAADCqR4cRt9utkpISud1u06Wccoy95429p45bYuw9cew9ddxSbI49JhawAgCA+NWjr4wAAADzCCMAAMAowggAADCKMAIAAIyK+zBSVlamzMxMJScnKycnRzt27Oiw/Ysvvqjhw4crOTlZI0eO1IYNG05Rpc4pLS3VJZdcon79+mnAgAGaNm2adu3a1WGf1atXy+VytXgkJyefooqdc99997Uax/DhwzvsEw9znpmZ2WrcLpdLc+fObbN9LM/3G2+8oSlTpigjI0Mul0vr1q1r8XvLslRcXKyBAweqT58+8vl8+uSTT056XrvvFSZ0NPbjx49rwYIFGjlypE477TRlZGSooKBA+/fv7/CckbxmTrWTzfnNN9/cagxXX331Sc8b63Muqc3Xvcvl0iOPPNLuOaNxzuM6jKxdu1aFhYUqKSlRdXW1srKylJeXp4MHD7bZfvv27Zo+fbpmzZqlnTt3atq0aZo2bZref//9U1x517z++uuaO3eu3n77bW3evFnHjx/XxIkT1dDQ0GG/lJQUHThwIPzYu3fvKarYWRdddFGLcbz55pvtto2XOf/nP//ZYsybN2+WJP3qV79qt0+szndDQ4OysrJUVlbW5u8ffvhh/fnPf1Z5ebneeecdnXbaacrLy9OxY8faPafd9wpTOhr70aNHVV1drUWLFqm6ulovv/yydu3apWuvvfak57XzmjHhZHMuSVdffXWLMTz//PMdnjMe5lxSizEfOHBAq1atksvl0i9+8YsOzxt1c27FsXHjxllz584N/9zc3GxlZGRYpaWlbba/7rrrrMmTJ7c4lpOTY/3mN7/p1jq728GDBy1J1uuvv95um6efftpKTU09dUV1k5KSEisrK6vT7eN1zufNm2ede+65VigUavP38TLfkqxXXnkl/HMoFLLS09OtRx55JHzsyJEjltvttp5//vl2z2P3vSIa/HDsbdmxY4clydq7d2+7bey+Zkxra9wzZsywpk6daus88TrnU6dOta688soO20TjnMftlZGmpiZVVVXJ5/OFjyUkJMjn86mysrLNPpWVlS3aS1JeXl677WNFXV2dJOnMM8/ssN0333yjIUOGyOv1aurUqfrggw9ORXmO++STT5SRkaFzzjlHN954o2pra9ttG49z3tTUpGeffVa//vWvO/xiyXiZ7/9tz5498vv9LeY0NTVVOTk57c5pJO8VsaKurk4ul0v9+/fvsJ2d10y02rZtmwYMGKDzzz9fc+bM0VdffdVu23id80AgoPXr12vWrFknbRttcx63YeTw4cNqbm6Wx+Npcdzj8cjv97fZx+/322ofC0KhkObPn68JEyZoxIgR7bY7//zztWrVKr366qt69tlnFQqFNH78eH3xxRensNquy8nJ0erVq7Vx40YtX75ce/bs0aWXXqr6+vo228fjnK9bt05HjhzRzTff3G6beJnvHzoxb3bmNJL3ilhw7NgxLViwQNOnT+/wy9Lsvmai0dVXX62//vWvqqio0EMPPaTXX39dkyZNUnNzc5vt43XOn3nmGfXr108///nPO2wXjXMeE9/ai8jNnTtX77///kk/D8zNzVVubm745/Hjx+uCCy7Qk08+qfvvv7+7y3TMpEmTwv88atQo5eTkaMiQIXrhhRc69X8L8WDlypWaNGmSMjIy2m0TL/ONth0/flzXXXedLMvS8uXLO2wbD6+Z66+/PvzPI0eO1KhRo3Tuuedq27ZtuuqqqwxWdmqtWrVKN95440kXo0fjnMftlZG0tDQlJiYqEAi0OB4IBJSent5mn/T0dFvto93tt9+u1157TVu3btXgwYNt9e3du7fGjBmj3bt3d1N1p0b//v01bNiwdscRb3O+d+9ebdmyRbfccoutfvEy3yfmzc6cRvJeEc1OBJG9e/dq8+bNtr9C/mSvmVhwzjnnKC0trd0xxNucS9I//vEP7dq1y/ZrX4qOOY/bMJKUlKTs7GxVVFSEj4VCIVVUVLT4P8L/LTc3t0V7Sdq8eXO77aOVZVm6/fbb9corr+jvf/+7hg4davsczc3Neu+99zRw4MBuqPDU+eabb/Tpp5+2O454mfMTnn76aQ0YMECTJ0+21S9e5nvo0KFKT09vMafBYFDvvPNOu3MayXtFtDoRRD755BNt2bJFP/rRj2yf42SvmVjwxRdf6Kuvvmp3DPE05yesXLlS2dnZysrKst03Kubc9Ara7rRmzRrL7XZbq1evtv79739bt956q9W/f3/L7/dblmVZN910k7Vw4cJw+7feesvq1auX9eijj1offvihVVJSYvXu3dt67733TA0hInPmzLFSU1Otbdu2WQcOHAg/jh49Gm7zw7EvXrzY2rRpk/Xpp59aVVVV1vXXX28lJydbH3zwgYkhROy3v/2ttW3bNmvPnj3WW2+9Zfl8PistLc06ePCgZVnxO+eW9d3dAGeffba1YMGCVr+Lp/mur6+3du7cae3cudOSZC1dutTauXNn+I6RBx980Orfv7/16quvWv/617+sqVOnWkOHDrW+/fbb8DmuvPJK67HHHgv/fLL3imjR0dibmpqsa6+91ho8eLBVU1PT4rXf2NgYPscPx36y10w06Gjc9fX11l133WVVVlZae/bssbZs2WJdfPHF1nnnnWcdO3YsfI54nPMT6urqrL59+1rLly9v8xyxMOdxHUYsy7Iee+wx6+yzz7aSkpKscePGWW+//Xb4d5dddpk1Y8aMFu1feOEFa9iwYVZSUpJ10UUXWevXrz/FFXedpDYfTz/9dLjND8c+f/788L8nj8djXXPNNVZ1dfWpL76L8vPzrYEDB1pJSUnWoEGDrPz8fGv37t3h38frnFuWZW3atMmSZO3atavV7+Jpvrdu3drmf98nxhcKhaxFixZZHo/Hcrvd1lVXXdXq38mQIUOskpKSFsc6eq+IFh2Nfc+ePe2+9rdu3Ro+xw/HfrLXTDToaNxHjx61Jk6caJ111llW7969rSFDhlizZ89uFSricc5PePLJJ60+ffpYR44cafMcsTDnLsuyrG699AIAANCBuF0zAgAAYgNhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFH/Hyv0TIIukHepAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(np.log10(hm.m), hm.satellite_occupation);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Density Profile\n", "\n", "The density profile from [arxiv:0906.3020](https://arxiv.org/abs/0906.3020) is already included as `PowerLawWithExpCut`:\n", "\n", "$$\n", "\\rho(r) = \\rho_s \\big(r/r_s \\big)^{-b}{\\rm exp}\\big[-a r/r_s\\big]\n", "$$\n", "\n", "and in this specific case we have $b=2$.\n", "\n", "However, the native way of defining density profile in `halomod` is to relate it to the characteristic scale $r_s$, which is related to the concentration parameter. Therefore, for each halo of different mass the shape of the density profile is different. But in this case we want to keep the shape of the profile the same for all halos. Although `halomod` does not provide a readily available solution, note:\n", "\n", "$$\n", "m \\sim r_s^3c^3(m,z)\n", "$$\n", "\n", "$$\n", "r_s \\sim m^{1/3}c^{-1}(m,z)\n", "$$\n", "\n", "Therefore, we only need to define a special concentration-mass relation to keep $r_s$ constant. Suppose we construct a C-M relation:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "from hmf.halos.mass_definitions import SOMean\n", "\n", "from halomod.concentration import CMRelation" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "class CMFlux(CMRelation):\n", " _defaults = {\"c_0\": 4}\n", " native_mdefs = (SOMean(),)\n", "\n", " def cm(self, m, z):\n", " return self.params[\"c_0\"] * (m * 10 ** (-11)) ** (1 / 3)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "hm = TracerHaloModel(\n", " halo_concentration_model=CMFlux,\n", " halo_profile_model=\"PowerLawWithExpCut\",\n", " halo_profile_params={\"b\": 2.0, \"a\": 1.0},\n", " hod_model=\"Constant\",\n", " transfer_model=\"EH\",\n", ")" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABGBUlEQVR4nO3deVxVdcI/8M+5K1xlEUEQRcANMxUUFbHGJSk0x9Q2pk1jytJHezSmabJFZ3pKe3oa6/dMPFlOatbMuJRZjWUpuSZqguC+oCxuoKgssl245/v7A71JAnKBy/cun/dr7gRfzrnnc7zS/XTuOd+jCCEEiIiIiCTRyA5ARERE7o1lhIiIiKRiGSEiIiKpWEaIiIhIKpYRIiIikoplhIiIiKRiGSEiIiKpWEaIiIhIKp3sAE2hqirOnTsHLy8vKIoiOw4RERE1gRACpaWlCA4OhkbT8PEPpygj586dQ0hIiOwYRERE1AynT59G165dG/y5U5QRLy8vALU74+3tLTkNERERNUVJSQlCQkKs7+MNcYoycv2jGW9vb5YRIiIiJ3OrUyx4AisRERFJxTJCREREUrGMEBERkVQsI0RERCQVywgRERFJxTJCREREUrGMEBERkVQsI0RERCQVywgRERFJZXMZ2bZtGyZMmIDg4GAoioJ169bdcp0tW7Zg0KBBMBqN6NmzJ5YvX96MqEREROSKbC4jZWVliIyMRHJycpOWz87Oxvjx4zF69GhkZGRgzpw5ePrpp/H999/bHJaIiIhcj833phk3bhzGjRvX5OUXL16M8PBw/PWvfwUA3HbbbdixYwfeffddxMfH27p5IiIicjF2v1Feamoq4uLi6ozFx8djzpw5Da5TVVWFqqoq6/clJSV2yfbGp0+gsCrfLs9NrkVB4zd5aqvtNZzi5p/UeQ6lgfEbvrZ+de2GVteXU65/Zb3RlXLDz278fwW1/9NcW7R2OY2i/LKOolwfvXbjrNqfK1CgKBrrPzWKBoqihcb6tQYajRYaRQeNokCr0UGjaKEoGmg1Wmg1umtjOuh1Bmg0Oui1Bui0euh1Buh1Ruh0Bhi0RhgMHjDoPeFp8ITBYILJox08je2g0+kb/NMlIvuyexnJz89HYGBgnbHAwECUlJSgoqICnp6eN62zcOFC/OUvf7F3NOyr2I/jRtXu2yFya6KBrx2MTgjoBKAXAjrg2tcK9EKBDgr0QgMdNNALLfSKDnroYFD0MChG6LVGeGg94aFrB099e5iM3jAZfeDTriM6eAXC37cLAjt2hU97P9m7SeSQ7F5GmmPu3LlISkqyfl9SUoKQkJBW305M+6HoUXm+1Z+XHIho23c/27fW9DVEQ8uK+pYR9f24weewjl/78xJ1fiLqXbbO/wsAiri2uqiz1o3f3fhPcW1FcW09ccOSdf6p/DKu3jCuQkBVrn8N688tSu3XFqV23AJAVYAaBbBAQY0C1NRzO/MapfZnlfUefxLXnskCoLreP0OIaz+qBlBe/yJ6IdBeFWinKminamESBrRTPNBO0x7t9T7w9eyEAO+u6OzXHWGdb0O3zr15xIbcgt3LSFBQEAoKCuqMFRQUwNvbu96jIgBgNBphNBrtHQ0vPrrE7tsgIsejWiww11Sh0lyO8soyVJkrUFlVhgpzOSqrrqKyqgxV1RWoNJehwlwGc005qqorUFVdDrOlElU1FaiyVKBarYLZUgWzqEKVMMMsqmFGDaoUCyoVFZUaFeUKUKZRIBQF1YqCK1oFV7QAoAKovPYoAnAGqDgEVAAoAHAE0AoBP4uAn0UHX5jQQdsBHT2CEOgTjl5dBiKy953waucr64+RqNXYvYzExsbi22+/rTO2ceNGxMbG2nvTRET10mi18NCa4GE0wdfL3+7bq6mpxoUr51B45SwuXDmLSyVnUXS1AMUVhSgxX0ZZTSnK1DJcFRUo0VSjWKuiSKuBRVFwUafgok4FcLX2UXMauPQzcGk1NJkCgTVAgMWIAE0HBJm6oVfnwYjtPx7BAaF23y+i1mJzGbl69SqysrKs32dnZyMjIwN+fn7o1q0b5s6di7Nnz2LFihUAgOnTp+P999/Hiy++iN///vf48ccfsXr1aqxfv7719oKIyIHpdHoEB4TaVBDKK8uQffYIcs4fxJmLx1BQmovLVRdQpJbgiqYC+ToV5RoNzuuB83ozgAKgpgA4/TNw+gMEVQt0sZgQbOiCnv4DMXrQwwjv0sd+O0nUAooQtn2ovmXLFowePfqm8alTp2L58uV48sknkZOTgy1bttRZ5/nnn8fhw4fRtWtXvPbaa3jyySebvM2SkhL4+PiguLgY3t7etsQlInJJqsWCrDMHkHliB3IuHkB+WS4uWApxVleBi7r6p5DqWi0QrnZET5/+uLPvJAzuOxoarbaNk5M7aer7t81lRAaWESKipss9dxy7Dn2L4/l7cbYiB3lKEU4bbj4xN6BGxW2WAAwMGomJv3kWAR2CJaQlV8YyQkREVrnnjiMlbSUOX0hFjuUcThksqL7hqiK9EIgwG9DXsy8mxvwHBvQeLjEtuQqWESIialBRaSG+2f537D27CUeUfJzX1z1y0rdKh6Hew/DImJd4Miw1G8sIERE1iWqxYNfBH5Cy/x84UHEQR4wW68+MqkCkuT3u7DIej93zIgwG+0+7QK6DZYSIiJpl//Gd+Dz1Pew1H65zrknnaoHRnkPwzPiF6OgbJDEhOQuWESIiahHVYsGmPavx3aFl2KM5ixJt7VU6XhYVdyIcT43+L0SED5SckhwZywgREbWaS0X5WLL+ZWyu2INz184v0QmBWLMfZo5ZhNt7DJackBwRywgREbU6s7kKKza8iQ35X+HYtRuNGlWBu0V3zJn4AQI7dpGckBwJywgREdnVV1s+wqfHk62lpEONit+a7sB/3v8ePIwmyenIETT1/bv+afqIiIhuYeKoZ7D6qXQ853sfgqsFrug0+NScigdWxGDlD+/KjkdOhEdGiIioxcrKS/HuFzPxXXWa9UTXu8yB+PPD/0IHnwDJ6UgWHhkhIqI2087khVefWIGV936DkebaOyH/aCjAI2tG45ttSyWnI0fHMkJERK0mJKg73p+2GX/q9Dg61qg4q1fw6qlFeGnpfSgrL5UdjxwUywgREbW6x8f9CZ+O/wbDqrygKgrWa7Px6D+GY9eBH2RHIwfEMkJERHYREtQdS57Zif/w/i28LSpOGYCkn5/H6o3/KzsaORiWESIisqsZkxfi76NWoFeVglKtBgvPfoS/rpwhOxY5EJYRIiKyu9u6R2PpIz8ipqo9ahQFy6t24A9/HwuzuUp2NHIALCNERNQmfL388dFTO/BbtQcA4Af9WTy97E5cvHJOcjKSjWWEiIjajEarxcLEdXi2fTwMqsA+j0o89Xk8Tp0+JDsaScQyQkREbW7WA+9gfo/n0cGiItsAzN7wO5w+f0J2LJKEZYSIiKS4b8RTeCv6HXSoUZFjAJ775gGcuZAjOxZJwDJCRETSDI8ch/+KWgAfi4qTRoHZ6yai4NJZ2bGojbGMEBGRVCOjJ+LPt/8ZXhYVx40qZn0xHpeK8mXHojbEMkJERNLFxTyEeX1eRnuLiqNGC2auHoei0kLZsaiNsIwQEZFDGDv8McztmQSTquKQsQbP/ZPzkLgLlhEiInIY9414Ci+FPQejKpDhUYVXPr1fdiRqAywjRETkUCaPno5En3EAgA26PHyw9k+SE5G9sYwQEZHDmXn//+Ce6i4AgKXF6/Hjns8lJyJ7YhkhIiKH9OYTa9G/Uo9KjYKF+//MSdFcGMsIERE5JA+jCf993yoEVQvk6xX86esEVFaVy45FdsAyQkREDiukcy/MHfBneKgCBzyq8eqnD8iORHbAMkJERA7trqEP4vc+4wEA3+vP4IMv50pORK2NZYSIiBzejPv/G/E1XQEAn135GsdzM+QGolbFMkJERE7h9cc+R3czUKLV4M0N06BaLLIjUSthGSEiIqdg8miH2ZHzoRMC6R6VSP7yj7IjUSthGSEiIqdx19AHca/oDQBYWfI9P65xESwjRETkVF555FN+XONiWEaIiMipmDzaYU7Un/lxjQthGSEiIqczesgDGC8iAPDjGlfAMkJERE7p5UdW8OMaF8EyQkRETunXH9f8/d/zZEeiZmIZISIipzV6yAOIV7sDAD6/+BXKykslJ6LmYBkhIiKn9ofJi+FXo+K8XsGiL2bIjkPNwDJCREROLaBDMO5rNwIA8F31PuSeOy45EdmKZYSIiJzec/e/hzAzUKrV4N1vZ8qOQzZiGSEiIqdnMBjxSLcnAQBbdeeRun+D3EBkE5YRIiJyCY/G/wFRlUbUKAo+TOWVNc6EZYSIiFzGs0P/DK0QSPOowJpN78uOQ03EMkJERC7jzoG/xYjqTgCAz7I/RE1NteRE1BQsI0RE5FLmjP0b2qkqThmA99f+QXYcagKWESIicindQ27HWE1/AMBXpSkoLSuSG4huiWWEiIhczvMP/B8CalQU6jT44OsXZcehW2AZISIil+PT3g93GaIBAJvKdqKyqlxyImoMywgREbmkGRPeho+ldpr4Jd+8KjsONYJlhIiIXFJH3yCM0vQBAHx/ZSOvrHFgzSojycnJCAsLg4eHB2JiYrBnz54Gl62ursbrr7+OHj16wMPDA5GRkdiwgTPjERGR/U2/968wqSpyDcCK7xbKjkMNsLmMrFq1CklJSZg/fz7S09MRGRmJ+Ph4XLhwod7lX331VXz44Yf429/+hsOHD2P69OmYPHky9u3b1+LwREREjenaKQwj1FAAwL/Pfw7VYpGciOqjCCGELSvExMRgyJAheP/92pntVFVFSEgInnvuObz00ks3LR8cHIxXXnkFM2f+cuOiBx54AJ6envjss8+atM2SkhL4+PiguLgY3t7etsQlIiI3dzx3Px798VFUaRS8GvwUEu6eIzuS22jq+7dNR0bMZjPS0tIQFxf3yxNoNIiLi0Nqamq961RVVcHDw6POmKenJ3bs2NHgdqqqqlBSUlLnQURE1By9QwdgeE3trKzrsldITkP1samMFBYWwmKxIDAwsM54YGAg8vPz610nPj4eixYtwokTJ6CqKjZu3Ii1a9fi/PnzDW5n4cKF8PHxsT5CQkJsiUlERFTH70e8Dp0QOGisxrc7WEgcjd2vpvl//+//oVevXujTpw8MBgNmzZqFxMREaDQNb3ru3LkoLi62Pk6fPm3vmERE5MKiIu5EjNkXALD6EG+g52hsKiP+/v7QarUoKCioM15QUICgoKB61wkICMC6detQVlaG3NxcHD16FO3bt0f37t0b3I7RaIS3t3edBxERUUs8PuQlKNfu6Ls9/WvZcegGNpURg8GA6OhopKSkWMdUVUVKSgpiY2MbXdfDwwNdunRBTU0NvvjiC0ycOLF5iYmIiJrhzoG/RXRVOwDAZ3vflpyGbmTzxzRJSUlYsmQJPvnkExw5cgQzZsxAWVkZEhMTAQBTpkzB3Llzrcvv3r0ba9euxalTp7B9+3aMHTsWqqrixRd5rwAiImpb99/2LABgr74IueeOS05D19lcRhISEvDOO+9g3rx5iIqKQkZGBjZs2GA9qTUvL6/OyamVlZV49dVX0bdvX0yePBldunTBjh074Ovr22o7QURE1BTj75iK7mbArFGwfNN82XHoGpvnGZGB84wQEVFreWfldHxS9RM6Vwt8++Q+6HR62ZFcll3mGSEiInJ2iWP/jHZq7Q30Vm16V3YcAssIERG5mY6+QYipqb0CdGPu55LTEMAyQkREbuj+Qf8JAMgwluPQyb2S0xDLCBERuZ2R0RNxW5UWFkXBZ9vekB3H7bGMEBGRWxrpfzcAIFU9gfLKMslp3BvLCBERuaUp8a/A16Likk6Dz75/S3Yct8YyQkREbsmrnS+GiTAAwJb89XLDuDmWESIicluP3FF7v5oDHtXYdeAH2XHcFssIERG5rUF9foMBVR4AgDW7/yo5jftiGSEiIrc2pkvtjVt3ac6gqLRQchr3xDJCRERu7bF7XkSnGhUlWg1WbOBlvjKwjBARkVszGIyI0fQGAOy6sk1yGvfEMkJERG5v0uCZAIBDBjOO52bIDeOGWEaIiMjtDe0fh95VGqiKgtXbFsmO43ZYRoiIiAAMbjcQALC3MkNuEDfEMkJERAQgYcQfoRUCJw0COzO/kx3HrbCMEBERAegecjv6XZtz5Ju0DySncS8sI0RERNfEdBwBAEhTT0G1WCSncR8sI0RERNc8GvciPFUV5/UK/v3Tctlx3AbLCBER0TUdfYMQVe0LANh07B9yw7gRlhEiIqIbjAi5DwCQri1AeWWZ5DTugWWEiIjoBg/e9Rw6WFQUazVYncI5R9oCywgREdENPIwmDFKDAQA7znwrOY17YBkhIiL6lXv6PAEAyDCUouDSWclpXB/LCBER0a+MjX0MwdUCVRoFK398W3Ycl8cyQkRE9CsarRaDtb0AAHsu75CcxvWxjBAREdVj4rU7+R40VuF47n7JaVwbywgREVE9hvaPQw+zAlVR8NXOZNlxXBrLCBERUQP66yMAAJkleyUncW0sI0RERA2Ij5wKADhsrEJ+4WnJaVwXywgREVED7hz4W3StFqhWFHyx7X9lx3FZLCNERESN6KeEAgDSL26XnMR1sYwQERE1YlTEQwCAg4ZSFF+9LDmNa2IZISIiakT8sMfgX6OiXKPBl1t5VY09sIwQERE1QqfTo78aCADYc3aj5DSuiWWEiIjoFoaH/hYAcEB3CWZzleQ0rodlhIiI6BbuG/EsvCwqirQafL3977LjuByWESIiolswebRD/5oOAICfsr+SnMb1sIwQERE1wZCguwAAB3AWqsUiOY1rYRkhIiJqggdGPgejKlCg12DTntWy47gUlhEiIqIm6OATgH5mEwBg85F/SU7jWlhGiIiImmigXywA4IAlW3IS18IyQkRE1ET3/2Y2tEIg1wDsObBJdhyXwTJCRETURCFB3XGbWQ8A+G4fL/FtLSwjRERENhjQLgoAsL/yiNwgLoRlhIiIyAbjhz4LAMgyWJB77rjkNK6BZYSIiMgGA3oNQ6gZUBUF/961RHYcl8AyQkREZKM+SlcAwP7CVMlJXAPLCBERkY2GhMYDAI7orqCmplpyGufHMkJERGSj397xFEyqiitazsbaGlhGiIiIbNTO5IXbzO0AADuOrZWcxvmxjBARETVDX+8oAMDR6iy5QVwAywgREVEzjB38ewDACYMFp8+fkJzGuTWrjCQnJyMsLAweHh6IiYnBnj17Gl3+vffeQ0REBDw9PRESEoLnn38elZWVzQpMRETkCAb0GoZu1y7x/SaVl/i2hM1lZNWqVUhKSsL8+fORnp6OyMhIxMfH48KFC/Uu/89//hMvvfQS5s+fjyNHjuDjjz/GqlWr8PLLL7c4PBERkUx9lGAAwP7CnZKTODeby8iiRYswbdo0JCYmom/fvli8eDFMJhOWLl1a7/I7d+7EHXfcgUcffRRhYWG455578Mgjj9zyaAoREZGjG9Kt9hLfw1pe4tsSNpURs9mMtLQ0xMXF/fIEGg3i4uKQmlr/xC/Dhw9HWlqatXycOnUK3377Le69994Gt1NVVYWSkpI6DyIiIkcz/o7fw1NVcUWnQcrPa2THcVo2lZHCwkJYLBYEBgbWGQ8MDER+fn696zz66KN4/fXXceedd0Kv16NHjx4YNWpUox/TLFy4ED4+PtZHSEiILTGJiIjahFc7X/QxmwAA2499ITmN87L71TRbtmzBggUL8H//939IT0/H2rVrsX79evzXf/1Xg+vMnTsXxcXF1sfp06ftHZOIiKhZbveKBAAcM5+UnMR56WxZ2N/fH1qtFgUFBXXGCwoKEBQUVO86r732Gp544gk8/fTTAID+/fujrKwMzzzzDF555RVoNDf3IaPRCKPRaEs0IiIiKeKjn8Rnu3bjhKEGZy7koGunMNmRnI5NR0YMBgOio6ORkpJiHVNVFSkpKYiNja13nfLy8psKh1arBQAIIWzNS0RE5FCiIu5EiFnAoij45qcPZcdxSjYdGQGApKQkTJ06FYMHD8bQoUPx3nvvoaysDImJiQCAKVOmoEuXLli4cCEAYMKECVi0aBEGDhyImJgYZGVl4bXXXsOECROspYSIiMiZ9VG64DTOIfPiDtlRnJLNZSQhIQEXL17EvHnzkJ+fj6ioKGzYsMF6UmteXl6dIyGvvvoqFEXBq6++irNnzyIgIAATJkzAm2++2Xp7QUREJNHgkHuwMX85jmovQ7VYoOF/bNtEEU7wWUlJSQl8fHxQXFwMb29v2XGIiIjqKC0rwpjVd6BCo8GiiFdx97AE2ZEcQlPfv3lvGiIiohbyaueLPtXXL/H9XHIa58MyQkRE1Apua9cPAHC8infxtRXLCBERUSv4ze0PAQBO6KtRVFooOY1zYRkhIiJqBcP7x6NjjQqzRsGG1E9lx3EqLCNEREStQKPVorfFFwCQfial8YWpDpYRIiKiVtLHZyAAIKsmT3IS58IyQkRE1EriBj0OADhpUHHuYq7kNM6DZYSIiKiVDOg1DMHVAqqiYMOu5bLjOA2WESIiolbUS/gDAPbnb5ecxHmwjBAREbWi2/1rbxybhYJbLEnXsYwQERG1orFDpwIAcg3A8dz9ktM4B5YRIiKiVhTepQ/CzbVf/7B3hdwwToJlhIiIqJX1VIIBAIcv7ZGcxDmwjBAREbWyyM4jAADHNZckJ3EOLCNERESt7N7hidAJgQK9BumHt8qO4/BYRoiIiFpZQIdg9DBrAQA/7v+X5DSOj2WEiIjIDnrpwwAAR0sy5QZxAiwjREREdjAo5G4AwHFdCVSLRXIax8YyQkREZAdjYx+HURW4otVga/o62XEcGssIERGRHXi180XvagMA4Kej6+SGcXAsI0RERHbS29gbAHC84qjkJI6NZYSIiMhOYnqOBwAc15ejsqpcchrHxTJCRERkJ2OGPAwvi4oyjQY/7OYlvg1hGSEiIrITg8GIXtUmAEB69kbJaRwXywgREZEd9TDVnjeSVZklOYnjYhkhIiKyo6E97gUAZBkqYDZXSU7jmFhGiIiI7GhU9GSY1NrzRrakrZUdxyGxjBAREdmRh9GEnmYPAMCerO8kp3FMLCNERER21sOjBwAgq5zzjdSHZYSIiMjOokPjAABZ+jLep6YeLCNERER2Fjf0ERhVgWKtBj9lrpcdx+GwjBAREdlZO5MXelbrAQA7j34tOY3jYRkhIiJqA931YQCAE1cPyQ3igFhGiIiI2kBk15EAgCxtCc8b+RWWESIiojZwz9DHoBMCl3QapB3dIjuOQ2EZISIiagMdfALQw6wFAGw/yMnPbsQyQkRE1Ea6a7sCAI4XZ0pO4lhYRoiIiNpI/86/AQCc1BTJDeJgWEaIiIjayD0xT0AjBPL1Cvaf2CU7jsNgGSEiImojgR27ILy69q13a+ZqyWkcB8sIERFRG+qudAYAHL2cLjmJ42AZISIiakP9AocDAE5pCiUncRwsI0RERG3o7sGPAgDO6BVk5R2UnMYxsIwQERG1oZDOvRBqrv06Je0fcsM4CJYRIiKiNtYDnQAAhwt/lpzEMbCMEBERtbE+/kMBAKdQIDmJY2AZISIiamNjBiUAAHL1Amcu5MgN4wBYRoiIiNpY79AodK4WEIqCzXtXyo4jHcsIERGRBOGqLwDgUP5OuUEcAMsIERGRBD28+gIAcmpOS04iH8sIERGRBDG9xwMATuqrUV5ZJjmNXCwjREREEtwReS+8LCoqNQq2pq2VHUcqlhEiIiIJdDo9uld7AgDSTm2UnEauZpWR5ORkhIWFwcPDAzExMdizZ0+Dy44aNQqKotz0GD9+fLNDExERuYJwYzgA4FTFMclJ5LK5jKxatQpJSUmYP38+0tPTERkZifj4eFy4cKHe5deuXYvz589bHwcPHoRWq8VDDz3U4vBERETObEDXkQCAU9qrUC0WyWnksbmMLFq0CNOmTUNiYiL69u2LxYsXw2QyYenSpfUu7+fnh6CgIOtj48aNMJlMLCNEROT24oYkQCcELuk0yDi2Q3YcaWwqI2azGWlpaYiLi/vlCTQaxMXFITU1tUnP8fHHH+N3v/sd2rVr1+AyVVVVKCkpqfMgIiJyNR18AhBurn0r3nHoS8lp5LGpjBQWFsJisSAwMLDOeGBgIPLz82+5/p49e3Dw4EE8/fTTjS63cOFC+Pj4WB8hISG2xCQiInIaoZogAMDxokzJSeRp06tpPv74Y/Tv3x9Dhw5tdLm5c+eiuLjY+jh9mhPCEBGRa7otYAgAIAeFkpPIY1MZ8ff3h1arRUFB3bsMFhQUICgoqNF1y8rKsHLlSjz11FO33I7RaIS3t3edBxERkSsaFVV707w8N75pnk1lxGAwIDo6GikpKdYxVVWRkpKC2NjYRtdds2YNqqqq8PjjjzcvKRERkQvqHTrA7W+aZ/PHNElJSViyZAk++eQTHDlyBDNmzEBZWRkSExMBAFOmTMHcuXNvWu/jjz/GpEmT0LFjx5anJiIiciHuftM8na0rJCQk4OLFi5g3bx7y8/MRFRWFDRs2WE9qzcvLg0ZTt+McO3YMO3bswA8//NA6qYmIiFxID6++2GlOddub5ilCCCE7xK2UlJTAx8cHxcXFPH+EiIhczta0rzDr4KvwUAW2PrIbJo+Gp79wJk19/+a9aYiIiCRz95vmsYwQERFJ5u43zWMZISIicgDWm+aVH5ecpO2xjBARETkA603zdKVud9M8lhEiIiIHUOemeSfc6xJflhEiIiIHcONN83464F4nsbKMEBEROYhfbpqXITdIG2MZISIichB9/AcDcL+b5rGMEBEROYgRAx4AUHvTvIJLZyWnaTssI0RERA7itu7R6FSjQlUUbE5bLTtOm2EZISIiciBhltpp0w+e2yE5SdthGSEiInIg4Z69AAC5VbmSk7QdlhEiIiIHMjD8LgBAtr7CbSY/YxkhIiJyIKMGPQCjKlCs1WD3Ife4Tw3LCBERkQNpZ/JCWLUWALD76HrJadoGywgREZGDCdV0BgBkFR+QnKRtsIwQERE5mIiAaABALi5JTtI2WEaIiIgczMgBDwJwn8nPWEaIiIgcTET4QARWX5/8bJXsOHbHMkJEROSAwtRrk5+d/UlyEvtjGSEiInJA4abeAIAcs+tPfsYyQkRE5IAGho0B4B6Tn7GMEBEROaCRgybDqAqUaDXYdfAH2XHsimWEiIjIAbUzeSG8WgcA2OXik5+xjBARETmo65OfnSw+KDmJfbGMEBEROajeAYMBALmKa09+xjJCRETkoEZFPQCgdvKz/MLTktPYD8sIERGRg+odGoXAahVCUfBj2mrZceyGZYSIiMiBhVlqJz87dM51Jz9jGSEiInJg4e0iAAC5Ljz5GcsIERGRAxsUXjv52Sl9JWpqqiWnsQ+WESIiIgc2Mvp+GFWBUq0Guw5ulB3HLlhGiIiIHJjJo5118rM9x1xz8jOWESIiIgcXqr0++dkhyUnsg2WEiIjIwfX2r538LM9FJz9jGSEiInJwv+k/CUDt5GcXr5yTG8YOWEaIiIgc3G3doxFQo0JVFGze+7nsOK2OZYSIiMgJhNV4AQAOntshOUnrYxkhIiJyAmGmHgCA3MpsyUlaH8sIERGRE4gMGQkAyNaVQ7VYJKdpXSwjRERETmBk9P3QCYErOg0yTuyUHadVsYwQERE5AV8vf4SZa9+2Uw9/LTlN62IZISIichLdlE4AgBOXM+QGaWUsI0RERE6iZ4cBAIA8tUByktbFMkJEROQkht02HgCQY1BRWlYkN0wrYhkhIiJyEtF9RsHHoqJaUbAlba3sOK2GZYSIiMhJaLRahFd7AgAyc7fIDdOKWEaIiIicSKgxFACQU3FCcpLWwzJCRETkRPoF3wEAyNGUSE7SelhGiIiInMioQQ9BEQIFeg2O52bIjtMqWEaIiIicSJB/CLpVKwCAbZlfSk7TOlhGiIiInEw34QcAOHZxr+QkrYNlhIiIyMn08O4LAMiznJOcpHU0q4wkJycjLCwMHh4eiImJwZ49expdvqioCDNnzkTnzp1hNBrRu3dvfPvtt80KTERE5O4G9xoLAMjWm2E2V0lO03I2l5FVq1YhKSkJ8+fPR3p6OiIjIxEfH48LFy7Uu7zZbMbdd9+NnJwcfP755zh27BiWLFmCLl26tDg8ERGRO4rtPxYmVUWFRoMdGd/IjtNiNpeRRYsWYdq0aUhMTETfvn2xePFimEwmLF26tN7lly5disuXL2PdunW44447EBYWhpEjRyIyMrLF4YmIiNyRwWBE92oDAGDvyR8kp2k5m8qI2WxGWloa4uLifnkCjQZxcXFITU2td52vv/4asbGxmDlzJgIDA9GvXz8sWLAAFoulZcmJiIjcWDddVwDAqdIjkpO0nE1lpLCwEBaLBYGBgXXGAwMDkZ+fX+86p06dwueffw6LxYJvv/0Wr732Gv7617/ijTfeaHA7VVVVKCkpqfMgIiKiX/TpNBQAkKdckZyk5ex+NY2qqujUqRM++ugjREdHIyEhAa+88goWL17c4DoLFy6Ej4+P9RESEmLvmERERE5lZOT9AIDTBgVnLuTIDdNCNpURf39/aLVaFBQU1BkvKChAUFBQvet07twZvXv3hlartY7ddtttyM/Ph9lsrneduXPnori42Po4ffq0LTGJiIhcXveQ29G5WgAAtqavlpymZWwqIwaDAdHR0UhJSbGOqaqKlJQUxMbG1rvOHXfcgaysLKiqah07fvw4OnfuDIPBUO86RqMR3t7edR5ERERUV6ha+/54+PwuyUlaxuaPaZKSkrBkyRJ88sknOHLkCGbMmIGysjIkJiYCAKZMmYK5c+dal58xYwYuX76M2bNn4/jx41i/fj0WLFiAmTNntt5eEBERuaEwU28AQK45T3KSltHZukJCQgIuXryIefPmIT8/H1FRUdiwYYP1pNa8vDxoNL90nJCQEHz//fd4/vnnMWDAAHTp0gWzZ8/Gn/70p9bbCyIiIjc0MOwurDyZhhx9BVSLBZobTolwJooQQsgOcSslJSXw8fFBcXExP7IhIiK6pqy8FCNWxcKsUfDxoHcxtH/crVdqQ019/+a9aYiIiJxUO5MXwqprj4bsOvpvyWmaz+aPaRyVxWJBdXW17BhkA61WC51OB0VRZEchInJa3TSBOI7zyCraLztKs7lEGbl69SrOnDkDJ/jEiX7FZDI1emUVERE1rpdfFDaVnkceCmVHaTanLyMWiwVnzpyByWRCQEAA/yvbSQghYDabcfHiRWRnZ6NXr151TnwmIqKmGX77ffhg13fI1asoKi2Er5e/7Eg2c/oyUl1dDSEEAgIC4OnpKTsO2cDT0xN6vR65ubkwm83w8PCQHYmIyOkM6BmLDjtUXNFpsDVtLSaOekZ2JJu5zH+K8oiIc+LRECKiltFotQivMQEAMk9vlZymefhOQERE5ORCPcIBADnlJyUnaR6WESIiIid3e/BwAECurlRykuZhGSEiInJyowY9BI0QuKDT4MipNNlxbMYyQs02efJkdOjQAQ8++KB17PTp0xg1ahT69u2LAQMGYM2aNRITEhG5h8COXdCtuvbcye0H1skN0wwsI9Rss2fPxooVK+qM6XQ6vPfeezh8+DB++OEHzJkzB2VlZZISEhG5j27wAwCcKOSREXIjo0aNgpeXV52xzp07IyoqCgAQFBQEf39/XL58WUI6IiL30sP7dgBAruWc5CS2YxlxQ9u2bcOECRMQHBwMRVGwbt26epdLTk5GWFgYPDw8EBMTgz179ti0nbS0NFgsFoSEhLRCaiIiasyQXmMBANn6GlRWlUtOYxuWETdUVlaGyMhIJCcnN7jMqlWrkJSUhPnz5yM9PR2RkZGIj4/HhQsXmrSNy5cvY8qUKfjoo49aKzYRETUidsA4tFNVVGoU7MhYLzuOTVhGJMrJyYGiKPjiiy8wYsQIeHp6YsiQIcjLy8P27dsxbNgwmEwmjBkzBkVFRa223XHjxuGNN97A5MmTG1xm0aJFmDZtGhITE9G3b18sXrwYJpMJS5cuveXzV1VVYdKkSXjppZcwfPjwVstNREQN0+n0CK82AgDST/4gOY1tWEYkyszMBAB88MEHWLBgAXbu3ImCggI8/vjjeOutt/D+++9j8+bNyMzMxLJly25af8GCBWjfvn2jj7y8PJtzmc1mpKWlIS4uzjqm0WgQFxeH1NTURtcVQuDJJ5/EXXfdhSeeeMLmbRMRUfOF6roCAE5dPSo5iW2c/t40vyaEQEW1Rcq2PfVam6alz8jIgJ+fH1atWoWOHTsCAEaOHIkdO3bg0KFDMJlqp/cdMmQI8vPzb1p/+vTpePjhhxvdRnBwsA17UKuwsBAWiwWBgYF1xgMDA3H06C9/wePi4pCZmYmysjJ07doVa9asgcViwapVqzBgwADruSiffvop+vfvb3MOIiKyTUSnIVh/KRt5yhXZUWzicmWkotqCvvO+l7Ltw6/Hw2Ro+h9pZmYmJk+ebC0iAJCXl4eEhARrEbk+NnHixJvW9/Pzg5+fX8tCt8CmTZvqHVdVtY2TEBERAIyMvB+LflyN0wYFZy7koGunMNmRmoQf00iUkZGBmJiYOmOZmZkYNmyY9fvKykocO3YMkZGRN61vr49p/P39odVqUVBQUGe8oKAAQUFBNj8fERG1je4ht6NztQAAbE1fLTlN07nckRFPvRaHX4+Xtu2mKikpQU5ODgYOHGgdy87ORnFxcZ2xAwcOQAhR78cc9vqYxmAwIDo6GikpKZg0aRKA2qMdKSkpmDVrls3PR0REbSdU9cZ5lOLw+V2yozSZy5URRVFs+qhElszMTGi1WvTr1886dv0cktDQ0DpjPXr0QPv27W96juZ+THP16lVkZWVZv8/OzrZuu1u3bgCApKQkTJ06FYMHD8bQoUPx3nvvoaysDImJiTZvj4iI2k6YqTd2WdKQa86VHaXJHP9d20VlZmYiIiICHh4edcZuPCpyfay+j2haYu/evRg9erT1+6SkJADA1KlTsXz5cgBAQkICLl68iHnz5iE/Px9RUVHYsGHDTSe1EhGRYxkYdhdWnkxDtr4SqsUCjbbpR+1lUYQQQnaIWykpKYGPjw+Ki4vh7e1d52eVlZXIzs5GeHh4nTd2cg58/YiIWldZeSlGrIqFWaNgyaC/Ylj/e6Rlaez9+0Y8gZWIiMiFtDN5Iby69mjI7qPOMRMrywgREZGL6aapvfLxZNEByUmahmWEiIjIxfTqOAgAkIuLkpM0DcsIERGRi/lN/0kAgFy9wKWim2fwdjQsI0RERC6mX88Y+NeosCgKNqd9ITvOLbGMEBERuaCwmtr5qQ6e2S45ya2xjBAREbmgMM/uAICcylOSk9waywgREZELGhAyEgCQrSuDapFzN/umYhkhIiJyQaMHPwidELis02B/VqrsOI1iGSEiInJBvl7+CK2ufZvfeehryWkaxzJCRETkokIRAAA4cTlDbpBbYBkhIiJyUT19BwAA8tQCyUkaxzJCzTZ58mR06NABDz74oHWsqKgIgwcPRlRUFPr164clS5ZITEhE5N5i+04EAGQbLCgtK5IbphEsI9Rss2fPxooVK+qMeXl5Ydu2bcjIyMDu3buxYMECXLp0SVJCIiL3NqjPb+BrUVHt4JOfsYxQs40aNQpeXl51xrRaLUwmEwCgqqoKQggIIWTEIyJyexqtFuE1tf9O3p+7RW6YRrCMuKFt27ZhwoQJCA4OhqIoWLduXb3LJScnIywsDB4eHoiJicGePXua9PxFRUWIjIxE165d8cc//hH+/v6tmJ6IiGwRaggDAGRXZMkN0giWETdUVlaGyMhIJCcnN7jMqlWrkJSUhPnz5yM9PR2RkZGIj4/HhQsXbvn8vr6+yMzMRHZ2Nv75z3+ioMCxT5wiInJltwcPBwDkakokJ2kYy4hEOTk5UBQFX3zxBUaMGAFPT08MGTIEeXl52L59O4YNGwaTyYQxY8agqKio1bY7btw4vPHGG5g8eXKDyyxatAjTpk1DYmIi+vbti8WLF8NkMmHp0qVN3k5gYCAiIyOxfbvj3xeBiMhVjY5+GBohUKDX4Fj2Ptlx6uV6ZUQIwFwm52HjuRGZmZkAgA8++AALFizAzp07UVBQgMcffxxvvfUW3n//fWzevBmZmZlYtmzZTesvWLAA7du3b/SRl5dn8x+h2WxGWloa4uLirGMajQZxcXFITW18Fr+CggKUlpYCAIqLi7Ft2zZERETYnIGIiFpHYMcu6FatAAC27l8rOU39dLIDtLrqcmBBsJxtv3wOMLRr8uIZGRnw8/PDqlWr0LFjRwDAyJEjsWPHDhw6dMh6IuiQIUOQn59/0/rTp0/Hww8/3Og2goNt/7MoLCyExWJBYGBgnfHAwEAcPXrU+n1cXBwyMzNRVlaGrl27Ys2aNdBqtXjmmWesJ64+99xz6N+/v80ZiIio9XQTHZGDSzheuFd2lHq5XhlxIpmZmZg8ebK1iABAXl4eEhISrEXk+tjEiRNvWt/Pzw9+fn5tkrU+mzZtqnc8IyOjbYMQEVGjevjcjm2V25BnOS87Sr1cr4zoTbVHKGRt2wYZGRmYO3dunbHMzEw8//zz1u8rKytx7NgxREZG3rT+ggULsGDBgka3cfjwYXTr1s2mXP7+/tBqtTedeFpQUICgoCCbnouIiOQb2nsclu3fhlP6GlRWlcPDaNv7lb25XhlRFJs+KpGlpKQEOTk5GDhwoHUsOzsbxcXFdcYOHDgAIUS9H3XY62Mag8GA6OhopKSkYNKkSQAAVVWRkpKCWbNm2fx8REQk17D+8fDa9yeUajXYlv4V7ol9RHakOlyvjDiJzMxMaLVa9OvXzzp2/RyS0NDQOmM9evRA+/btb3qO5n5Mc/XqVWRl/XK9eXZ2tnXb14+iJCUlYerUqRg8eDCGDh2K9957D2VlZUhMTLR5e0REJJdOp0d4tQf2a81IP7WRZYRqZWZmIiIiAh4eHnXGbjwqcn2svo9oWmLv3r0YPXq09fukpCQAwNSpU7F8+XIAQEJCAi5evIh58+YhPz8fUVFR2LBhw00ntRIRkXPoZgjBfpzEqbJjsqPcRBFOMFd3SUkJfHx8UFxcDG9v7zo/q6ysRHZ2NsLDw+u8sZNz4OtHRNQ2/rHhbbxV8Ck6Vwv88PTBNtlmY+/fN3K9eUaIiIjoJqOjE6AIgfN6BVl5bVNGmoplhIiIyA0EB4RaJz/bkrFGcpq6WEaIiIjcRChqL3o4erFpNz5tKywjREREbqKHd+0VnLkONvkZywgREZGbGNbntwCAbH0NyspLJaf5RbPKSHJyMsLCwuDh4YGYmBjs2dPw4Z7ly5dDUZQ6D141QURE1PaG9bsH3hYVVRoFW9O/lB3HyuYysmrVKiQlJWH+/PlIT09HZGQk4uPjceHChQbX8fb2xvnz562P3NzcFoUmIiIi22m0WnSv9gQApGfXf38xGWwuI4sWLcK0adOQmJiIvn37YvHixTCZTFi6dGmD6yiKgqCgIOuDE2cRERHJEWqsneU7u+KE5CS/sKmMmM1mpKWlIS4u7pcn0GgQFxeH1NTUBte7evUqQkNDERISgokTJ+LQoUONbqeqqgolJSV1HkRERNRy/bv8BgCQrXWc91abykhhYSEsFstNRzYCAwORn59f7zoRERFYunQpvvrqK3z22WdQVRXDhw/HmTNnGtzOwoUL4ePjY32EhITYEpOIiIgacNfgh6ERAhd1GhzM2i07DoA2uJomNjYWU6ZMQVRUFEaOHIm1a9ciICAAH374YYPrzJ07F8XFxdbH6dOn7R2TiIjILQR0CEbotcnPth9YJzfMNTaVEX9/f2i1WhQUFNQZLygoQFBQUJOeQ6/XY+DAgXXuGvtrRqMR3t7edR7keCZPnowOHTrgwQcfvOln5eXlCA0NxQsvvCAhGRERNSYMAQCA45fSJCepZVMZMRgMiI6ORkpKinVMVVWkpKQgNja2Sc9hsVhw4MABdO7c2bak5HBmz56NFStW1PuzN998E8OGDWvjRERE1BQ9fQcAAHLV+k+xaGs2f0yTlJSEJUuW4JNPPsGRI0cwY8YMlJWVITExEQAwZcoUzJ0717r866+/jh9++AGnTp1Ceno6Hn/8ceTm5uLpp59uvb0gKUaNGgUvL6+bxk+cOIGjR49i3LhxElIREdGtDO87CQCQbVBRfPWy3DBoRhlJSEjAO++8g3nz5iEqKgoZGRnYsGGD9aTWvLw8nD//yzSzV65cwbRp03Dbbbfh3nvvRUlJCXbu3Im+ffu23l6QTbZt24YJEyYgODgYiqJg3bp19S5ny+R2N3rhhRewcOHCVkxMREStaVCf36CDRUWNomDL3s9lx2neCayzZs1Cbm4uqqqqsHv3bsTExFh/tmXLFixfvtz6/bvvvmtdNj8/H+vXr8fAgQNbHJyar6ysDJGRkUhOTm5wmeZMbgcAX331FXr37o3evXu3dmwiImoltZOftQMAZORtlpyG96aRKicnB4qi4IsvvsCIESPg6emJIUOGIC8vD9u3b8ewYcNgMpkwZswYFBUVtdp2x40bhzfeeAOTJ09ucJnmTG4HALt27cLKlSsRFhaGF154AUuWLMHrr7/eatmJiKh1hHl0BwDkVJySnIRlRKrMzEwAwAcffIAFCxZg586dKCgowOOPP4633noL77//PjZv3ozMzEwsW7bspvUXLFiA9u3bN/rIy8uzOVdzJ7cDaueIOX36NHJycvDOO+9g2rRpmDdvns0ZiIjIvvqHjAAAnNJdhWqxSM2ik7p1OxBCoKKmQsq2PXWeUBSlyctnZGTAz88Pq1atQseOHQEAI0eOxI4dO3Do0CGYTCYAwJAhQ+qdVG769Ol4+OGHG91GcHCwDXtQq7HJ7Y4ePWr9Pi4uDpmZmSgrK0PXrl2xZs2aJl9VRUREco0Z8jDeyPs/XNZpkHFiJwb1+Y20LC5XRipqKhDzz5hbL2gHux/dDZPe1OTlMzMzMXnyZGsRAWpPAE5ISLAWketjEydOvGl9Pz8/+Pn5tSx0C2za1PhNlp588sm2CUJERDbz9fJHmFmDLKPATwe/lFpG+DGNRBkZGXVO/gVqC8qN83NUVlbi2LFjiIyMvGl9e31M0xqT2xERkeML1dQeAc+6kik1h8sdGfHUeWL3o3Lm2vfUeTZ52ZKSEuTk5NS5sig7OxvFxcV1xg4cOAAhBPr373/Tc9jrY5obJ7ebNGkSgF8mt5s1a5bNz0dERI6pt99ApJR+hxw0fqWkvblcGVEUxaaPSmTJzMyEVqtFv379rGPXzyEJDQ2tM9ajRw+0b9/+pudo7sc0V69erTMdf3Z2tnXb3bp1A1A7ud3UqVMxePBgDB06FO+9916dye2IiMj53dn/fnyw8zvk6gUuFeWjo6+co98uV0acRWZmJiIiIuDh4VFn7NdzsGRmZtb7EU1L7N27F6NHj7Z+n5SUBACYOnWqdY6YhIQEXLx4EfPmzUN+fj6ioqLqTG5HRETOb0CvYfDfpqJQp0HKz6vx8N3/KSWHIoQQUrZsg5KSEvj4+KC4uPimm+ZVVlYiOzsb4eHhdd7YyTnw9SMikuvPn/wOKiyYFPOfrX4Sa2Pv3zfikREiIiI39uepK2VH4NU0REREJBfLCBEREUnFMkJERERSsYwQERGRVCwjREREJJXLlBEnuEKZ6qGqquwIREQkmdNf2qvX66EoCi5evIiAgACb7ppL8gghYDabcfHiRWg0GhgMBtmRiIhIEqcvI1qtFl27dsWZM2eQk5MjOw7ZyGQyoVu3btBoXOYgHRER2cjpywgAtG/fHr169UJ1dbXsKGQDrVYLnU7Ho1lERG7OJcoIUPvGptVqZccgIiIiG/HYOBEREUnFMkJERERSsYwQERGRVE5xzsj1OURKSkokJyEiIqKmuv6+fau5wJyijJSWlgIAQkJCJCchIiIiW5WWlsLHx6fBnyvCCaYuVVUV586dg5eXV6teBlpSUoKQkBCcPn0a3t7erfa8joT76BrcYR8B99hP7qNr4D42jRACpaWlCA4ObnQ+Kac4MqLRaNC1a1e7Pb+3t7fL/mW6jvvoGtxhHwH32E/uo2vgPt5aY0dEruMJrERERCQVywgRERFJ5dZlxGg0Yv78+TAajbKj2A330TW4wz4C7rGf3EfXwH1sXU5xAisRERG5Lrc+MkJERETysYwQERGRVCwjREREJBXLCBEREUnlVmUkJycHTz31FMLDw+Hp6YkePXpg/vz5MJvNja5XWVmJmTNnomPHjmjfvj0eeOABFBQUtFFq27355psYPnw4TCYTfH19m7TOk08+CUVR6jzGjh1r36At0Jx9FEJg3rx56Ny5Mzw9PREXF4cTJ07YN2gLXL58GY899hi8vb3h6+uLp556ClevXm10nVGjRt30Ok6fPr2NEjdNcnIywsLC4OHhgZiYGOzZs6fR5desWYM+ffrAw8MD/fv3x7fffttGSZvPln1cvnz5Ta+Zh4dHG6a1zbZt2zBhwgQEBwdDURSsW7fuluts2bIFgwYNgtFoRM+ePbF8+XK752wJW/dxy5YtN72GiqIgPz+/bQI3w8KFCzFkyBB4eXmhU6dOmDRpEo4dO3bL9ez1++hWZeTo0aNQVRUffvghDh06hHfffReLFy/Gyy+/3Oh6zz//PL755husWbMGW7duxblz53D//fe3UWrbmc1mPPTQQ5gxY4ZN640dOxbnz5+3Pv71r3/ZKWHLNWcf3377bfzv//4vFi9ejN27d6Ndu3aIj49HZWWlHZM232OPPYZDhw5h48aN+Pe//41t27bhmWeeueV606ZNq/M6vv32222QtmlWrVqFpKQkzJ8/H+np6YiMjER8fDwuXLhQ7/I7d+7EI488gqeeegr79u3DpEmTMGnSJBw8eLCNkzedrfsI1M5weeNrlpub24aJbVNWVobIyEgkJyc3afns7GyMHz8eo0ePRkZGBubMmYOnn34a33//vZ2TNp+t+3jdsWPH6ryOnTp1slPCltu6dStmzpyJXbt2YePGjaiursY999yDsrKyBtex6++jcHNvv/22CA8Pb/DnRUVFQq/XizVr1ljHjhw5IgCI1NTUtojYbMuWLRM+Pj5NWnbq1Kli4sSJds1jD03dR1VVRVBQkPif//kf61hRUZEwGo3iX//6lx0TNs/hw4cFAPHzzz9bx7777juhKIo4e/Zsg+uNHDlSzJ49uw0SNs/QoUPFzJkzrd9bLBYRHBwsFi5cWO/yDz/8sBg/fnydsZiYGPHss8/aNWdL2LqPtvyeOhoA4ssvv2x0mRdffFHcfvvtdcYSEhJEfHy8HZO1nqbs4+bNmwUAceXKlTbJZA8XLlwQAMTWrVsbXMaev49udWSkPsXFxfDz82vw52lpaaiurkZcXJx1rE+fPujWrRtSU1PbImKb2bJlCzp16oSIiAjMmDEDly5dkh2p1WRnZyM/P7/O6+jj44OYmBiHfB1TU1Ph6+uLwYMHW8fi4uKg0Wiwe/fuRtf9xz/+AX9/f/Tr1w9z585FeXm5veM2idlsRlpaWp3XQKPRIC4ursHXIDU1tc7yABAfH++QrxnQvH0EgKtXryI0NBQhISGYOHEiDh061BZx24SzvYYtERUVhc6dO+Puu+/GTz/9JDuOTYqLiwGg0fdDe76WTnGjPHvJysrC3/72N7zzzjsNLpOfnw+DwXDTeQmBgYEO/XmgrcaOHYv7778f4eHhOHnyJF5++WWMGzcOqamp0Gq1suO12PXXKjAwsM64o76O+fn5Nx3i1el08PPzazTvo48+itDQUAQHB2P//v3405/+hGPHjmHt2rX2jnxLhYWFsFgs9b4GR48erXed/Px8p3nNgObtY0REBJYuXYoBAwaguLgY77zzDoYPH45Dhw7Z9QahbaWh17CkpAQVFRXw9PSUlKz1dO7cGYsXL8bgwYNRVVWFv//97xg1ahR2796NQYMGyY53S6qqYs6cObjjjjvQr1+/Bpez5++jSxwZeemll+o9eejGx6//RXD27FmMHTsWDz30EKZNmyYpedM1Zx9t8bvf/Q733Xcf+vfvj0mTJuHf//43fv75Z2zZsqX1duIW7L2PjsDe+/jMM88gPj4e/fv3x2OPPYYVK1bgyy+/xMmTJ1txL6g1xcbGYsqUKYiKisLIkSOxdu1aBAQE4MMPP5QdjZooIiICzz77LKKjozF8+HAsXboUw4cPx7vvvis7WpPMnDkTBw8exMqVK6VlcIkjI3/4wx/w5JNPNrpM9+7drV+fO3cOo0ePxvDhw/HRRx81ul5QUBDMZjOKiorqHB0pKChAUFBQS2LbxNZ9bKnu3bvD398fWVlZGDNmTKs9b2PsuY/XX6uCggJ07tzZOl5QUICoqKhmPWdzNHUfg4KCbjrhsaamBpcvX7bp711MTAyA2qOAPXr0sDlva/L394dWq73pSrTGfpeCgoJsWl625uzjr+n1egwcOBBZWVn2iNjmGnoNvb29XeKoSEOGDh2KHTt2yI5xS7NmzbKeIH+rI3H2/H10iTISEBCAgICAJi179uxZjB49GtHR0Vi2bBk0msYPDkVHR0Ov1yMlJQUPPPAAgNozpvPy8hAbG9vi7E1lyz62hjNnzuDSpUt13rjtzZ77GB4ejqCgIKSkpFjLR0lJCXbv3m3zVUct0dR9jI2NRVFREdLS0hAdHQ0A+PHHH6GqqrVgNEVGRgYAtOnr2BCDwYDo6GikpKRg0qRJAGoPD6ekpGDWrFn1rhMbG4uUlBTMmTPHOrZx48Y2/d2zRXP28dcsFgsOHDiAe++9145J205sbOxNl3868mvYWjIyMhzi964hQgg899xz+PLLL7FlyxaEh4ffch27/j62+BRYJ3LmzBnRs2dPMWbMGHHmzBlx/vx56+PGZSIiIsTu3butY9OnTxfdunUTP/74o9i7d6+IjY0VsbGxMnahSXJzc8W+ffvEX/7yF9G+fXuxb98+sW/fPlFaWmpdJiIiQqxdu1YIIURpaal44YUXRGpqqsjOzhabNm0SgwYNEr169RKVlZWydqNRtu6jEEK89dZbwtfXV3z11Vdi//79YuLEiSI8PFxUVFTI2IVbGjt2rBg4cKDYvXu32LFjh+jVq5d45JFHrD//9d/VrKws8frrr4u9e/eK7Oxs8dVXX4nu3buLESNGyNqFm6xcuVIYjUaxfPlycfjwYfHMM88IX19fkZ+fL4QQ4oknnhAvvfSSdfmffvpJ6HQ68c4774gjR46I+fPnC71eLw4cOCBrF27J1n38y1/+Ir7//ntx8uRJkZaWJn73u98JDw8PcejQIVm70KjS0lLr7xsAsWjRIrFv3z6Rm5srhBDipZdeEk888YR1+VOnTgmTyST++Mc/iiNHjojk5GSh1WrFhg0bZO3CLdm6j++++65Yt26dOHHihDhw4ICYPXu20Gg0YtOmTbJ24ZZmzJghfHx8xJYtW+q8F5aXl1uXacvfR7cqI8uWLRMA6n1cl52dLQCIzZs3W8cqKirEf/zHf4gOHToIk8kkJk+eXKfAOJqpU6fWu4837hMAsWzZMiGEEOXl5eKee+4RAQEBQq/Xi9DQUDFt2jTrvzwdka37KETt5b2vvfaaCAwMFEajUYwZM0YcO3as7cM30aVLl8Qjjzwi2rdvL7y9vUViYmKdsvXrv6t5eXlixIgRws/PTxiNRtGzZ0/xxz/+URQXF0vag/r97W9/E926dRMGg0EMHTpU7Nq1y/qzkSNHiqlTp9ZZfvXq1aJ3797CYDCI22+/Xaxfv76NE9vOln2cM2eOddnAwEBx7733ivT0dAmpm+b6Zay/flzfp6lTp4qRI0fetE5UVJQwGAyie/fudX4vHZGt+/jf//3fokePHsLDw0P4+fmJUaNGiR9//FFO+CZq6L3wxtemLX8flWuhiIiIiKRwiatpiIiIyHmxjBAREZFULCNEREQkFcsIERERScUyQkRERFKxjBAREZFULCNEREQkFcsIERERScUyQkRERFKxjBAREZFULCNEREQkFcsIERERSfX/AWdYMXgAwopmAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(np.log10(hm.k_hm), hm.tracer_profile.u(hm.k_hm, m=1e12), label=\"$m = 10^{12}$\")\n", "plt.plot(np.log10(hm.k_hm), hm.tracer_profile.u(hm.k_hm, m=1e13), label=\"$m = 10^{13}$\")\n", "plt.plot(np.log10(hm.k_hm), hm.tracer_profile.u(hm.k_hm, m=1e14), label=\"$m = 10^{14}$\")\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One can see that indeed the density profile is now independant of halo mass" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tuning parameters" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So far the parameters are randomly set without clear physical meanings. We can easily tune these parameters to desired physical values.\n", "\n", "Suppose we want the mass density of galaxies to be $10^{-2}$ of the total critical density:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "rhoc = hm.cosmo.critical_density0.to(\"Msun/Mpc^3\").value * hm.cosmo.h**2" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "np.float64(6.046025739923064e-13)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hm.mean_tracer_den / rhoc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That means the parameter `logA` for the HOD should be changed to:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "np.float64(12.218530008219606)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "-np.log10(hm.mean_tracer_den / rhoc)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can simply set this on the existing model (everything that's dependent on it will be auto-updated):" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "hm.hod_params = {\"logA\": -np.log10(hm.mean_tracer_den / rhoc)}" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "np.float64(1.000000000000003)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hm.mean_tracer_den / rhoc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The density profile should satisfy $r_s/a = \\lambda_{\\rm mfp}$. $r_s$ can be obtained as:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "rs = hm.halo_profile.scale_radius(1e11)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Just to make sure, we calculate $r_s$ for a different halo mass:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "np.float64(0.027958379670006795)" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hm.halo_profile.scale_radius(1e12)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "in the units of Mpc/h. Assume we want $\\lambda_{\\rm mfp} = 10$ Mpc/h:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "hm.halo_profile_params = {\"a\": rs / 10}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check the density profile to see the cut-off:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG1CAYAAAAFuNXgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABS3klEQVR4nO3de1xTdeMH8M82xsZ1gMj9qngXRVEQ7yZlWpo9VmalZqbVY2VRPmVlZv1Su1hW+mSZmtaTWt4qrxneFUVBzQzxAshNboKM6xjb+f2BrUhEhhtnG5/367WX7Oyc7TNPtI9n33O+EkEQBBARERHZCKnYAYiIiIhMieWGiIiIbArLDREREdkUlhsiIiKyKSw3REREZFNYboiIiMimsNwQERGRTWG5ISIiIptiJ3aAlqbX65GbmwsXFxdIJBKx4xAREVETCIKAsrIy+Pn5QSpt/NhMqys3ubm5CAwMFDsGERERNUNWVhYCAgIaXafVlRsXFxcAdX85rq6uIqchIiKiplCr1QgMDDR8jjem1ZWbP7+KcnV1ZbkhIiKyMk0ZUsIBxURERGRTWG6IiIjIprDcEBERkU1huSEiIiKbwnJDRERENoXlhoiIiGwKyw0RERHZFJYbIiIisimilpsDBw5g9OjR8PPzg0QiwZYtW265zb59+9C7d28oFAqEhYXh66+/NntOIiIish6ilpuKigr07NkTS5cubdL66enpuOeeezBs2DCcOnUKL7zwAp588kns2rXLzEmJiIjIWog6/cLIkSMxcuTIJq+/bNkyhIaGYtGiRQCALl264NChQ/j4448xYsQIc8UkIiIiK2JVY24SEhIQGxtbb9mIESOQkJBw0200Gg3UanW9GxEREdkuq5o4My8vD97e3vWWeXt7Q61Wo6qqCg4ODjdss2DBAsybN8/s2dKLKvD65jNwtLeDs0IGJ4UdnBV2cLp+c1XawdVBDlelHK4OdnBVyqFylMNFYdekScCIiIioaayq3DTH7NmzERcXZ7j/55Tppna1XIMjl64avZ1cJoG7oz08nOpu7k72aOusQFsXBbxcFPB2VcLLVQFvFyXcHOUsQkRERLdgVeXGx8cH+fn59Zbl5+fD1dW1waM2AKBQKKBQKMyeLcTTCZ88HIEKjQ4VmlqUa2pRWVOLco0O5ZpalFVroa7SorRKC3V1LdRVWmhq9dDqBBSUaVBQprnlazjIZfB3d0CA4eaIAHcHtPN0RqinExzsZWZ/n0RERJbOqspNTEwMtm/fXm/Z7t27ERMTI1Kiv3g6K3BfhL9R21RrdSiuqDHcSiprcLW8BkXlGuSrNSgoq0ZhmQb56mqUVGpRpdXhYkE5LhaU3/BcEgngp3JAey9ntG/rhA5eLujm54pOPi5Qyll6iIio9RC13JSXl+PixYuG++np6Th16hQ8PDwQFBSE2bNnIycnB2vWrAEAPP3001iyZAn+85//4IknnsCePXvw/fffY9u2bWK9hduilMvg5+YAP7eGjzr9naZWhyvXqpFdUoXskkrDn5nFlbhUWIHSKi1yrlUh51oVDpwvNGwnk0rQwcsZXf1c0d1PhYggN4T7qyCXWdVYciIioiaTCIIgiPXi+/btw7Bhw25YPnnyZHz99dd4/PHHkZGRgX379tXb5sUXX8Qff/yBgIAAzJkzB48//niTX1OtVkOlUqG0tBSurq4meBfiEwQBxRU1SCuqwKWCclwqLMe5vDKczVWjuKLmhvWVcil6B7mjb4gHokI90CvIDY72VnUQj4iIWhljPr9FLTdisMVyczOCICBPXY2zOWr8nluK33NKkXS5BCWV2nrr2UkliAr1wB2dvTC8izdCPZ1ESkxERNQwlptGtKZy0xC9XsClwnIkZhQjMb0Yx9OLkVtaXW+ddp5OuKOzF2K7eiMqxANSKc/QIiIicbHcNKK1l5uGZBRVYM+5Auw5V4Bj6Veh1f31n4SfSomxvfzxr97+CPNyETElERG1Ziw3jWC5aVxZtRaHLhTh15QC7P4jD+rqWsNjPQJUuL+XP+6L8IeHk72IKYmIqLVhuWkEy03TVWt12HOuAJuSs7EvtRC1+rr/VBR2UvyrdwCmDgxFmJezyCmJiKg1YLlpBMtN81wt1+Dn07n4ISkbZ3P/mp9reGcvTB0Uiph2bXj1ZCIiMhuWm0aw3NweQRCQmF6M5QfTEX8uH3/+19Pd3xXP3dEBd3X1ZskhIiKTY7lpBMuN6aQVlmPFoXRsTM5GtVYPAOgb4o7Zo7qgd5C7yOmIiMiWsNw0guXG9IoravDVwTSsOJQOTW1dyRkV7oNZIzrzmjlERGQSLDeNYLkxnyulVfh493n8kJQNQai7OOBj/YLx4p0doXKQix2PiIisGMtNI1huzO9cnhoLd5zDvtS6Oa58XJVYMC4cwzp5iZyMiIislTGf35w9kUyus48rvp4Shf89GY2QNo7IU1djyqrjmPXDaZRWaW/9BERERLeB5YbMZkCYJ3bMHIypA0MhkQA/JGXjro/3Y8+5fLGjERGRDWO5IbNysJdhzr1dseHpGLTzdEK+WoMnvj6BWT+cRlWNTux4RERkg1huqEVEBntg+8xBmD64neEozrjPjyCruFLsaEREZGNYbqjFKOUyvDaqC9ZO6wdPZ3v8cUWN0UsO4eCFQrGjERGRDWG5oRbXr10b/PzcQPQMUOFapRaTVybii/2X0MpO3CMiIjNhuSFR+KocsP6pGDwYGQC9ACzYcQ7Prj2JypraW29MRETUCJYbEo1SLsP7D/TAO2O7w04qwbbfrmD8F0dRUlEjdjQiIrJiLDckKolEgon9grF2ej94ONnjTE4pxn+ZgAJ1tdjRiIjISrHckEXoG+KB9dP7wdtVgfP55XjoiwRkl/BMKiIiMh7LDVmMDt4u+OGp/ghwd0DG1Uo8tCwBaYXlYsciIiIrw3JDFiWojSM2PN0f7ds6Ibe0Gg99cRQpV9RixyIiIivCckMWx0elxPqnYtDV1xVF5Ro8/OVR/J5TKnYsIiKyEiw3ZJE8nRVYO70fege5obRKi8dXJeLy1QqxYxERkRVguSGLpXKQY/UTUdeP4NRg8spEFJVrxI5FREQWjuWGLJqLUo6vp/Q1DDKe+vVxVGh4oT8iIro5lhuyeF6uSqx+IgrujnKczi7FjO+SodXpxY5FREQWiuWGrEL7ts5Y8XhfKOVS7EstxOxNZzgXFRERNYjlhqxG7yB3LH2kN2RSCTYkZePDX1LFjkRERBaI5YasyvAu3nh3bHcAwNK9l/DjqRyRExERkaVhuSGr83BUEGYMaw8AeHXjGZzPLxM5ERERWRKWG7JKcXd2woCwNqjS6vD0t0ko5xlURER0HcsNWSWZVIJPH+4FH1cl0gor8MqG3zjAmIiIALDckBVr46zA0kd7w04qwbYzV7DqcIbYkYiIyAKw3JBViwx2x+v3dAEAzN+egqTLxSInIiIisbHckNV7vH8I7unhi1q9gH//L5lTNBARtXIsN2T1JBIJ3hvXA+3bOiFfrcEL605Br+f4GyKi1orlhmyCs8IOyx6LhFIuxaGLRfjfsctiRyIiIpGw3JDN6ODtglfv7gwAmL/9HDKvVoqciIiIxMByQzZlUkwI+rXzQJVWh5c3nObXU0RErRDLDdkUqVSCDx7oCUd7GRLTi/H1kQyxIxERUQtjuSGbE+jhiNdG1Z0e/v6uc0grLBc5ERERtSSWG7JJj0YHYWCYJ6q1erz8w2no+PUUEVGrwXJDNkkikeC9B3rARWGH5Mxr+OpgmtiRiIiohbDckM3yd3PAnHu7AgAW7T6PC5w9nIioVWC5IZv2YJ8ADOvUFjW1ery2+Qwn1yQiagVYbsimSSQSvHt/OBzkMhzPKMGWUzliRyIiIjNjuSGb5+fmgGfvCANQd3G/smqtyImIiMicWG6oVXhyUChCPZ1QWKbB4l8viB2HiIjMiOWGWgWFnQxvjekGAPj6SAZS8zi4mIjIVrHcUKsxpGNbjOjmDZ1ewJs//s7BxURENorlhlqVOfd2hVIuxbH0Yvx0OlfsOEREZAYsN9SqBLg7YsbQPwcXp6BcUytyIiIiMjWWG2p1pg1uh+A2jshXa/BZPAcXExHZGpYbanWUchnmjq67cvGKQ+m4WMCJNYmIbAnLDbVKd3T2xvDOXqjVC/hg1zmx4xARkQmx3FCr9erIzpBKgF1n85F0uUTsOEREZCKil5ulS5ciJCQESqUS0dHRSExMbHT9xYsXo1OnTnBwcEBgYCBefPFFVFdXt1BasiUdvF3wYGQgAOC9Hed4ajgRkY0QtdysX78ecXFxmDt3LpKTk9GzZ0+MGDECBQUFDa7/3Xff4dVXX8XcuXORkpKCFStWYP369XjttddaODnZihfu7ACFnRSJGcXYm9rwf3dERGRdRC03H330EaZNm4YpU6aga9euWLZsGRwdHbFy5coG1z9y5AgGDBiARx55BCEhIbjrrrswYcKEWx7tIboZX5UDHh8QAgB4b0cqdHoevSEisnailZuamhokJSUhNjb2rzBSKWJjY5GQkNDgNv3790dSUpKhzKSlpWH79u0YNWrUTV9Ho9FArVbXuxH93b+HhMFVaYfU/DJsPslZw4mIrJ1o5aaoqAg6nQ7e3t71lnt7eyMvL6/BbR555BG8/fbbGDhwIORyOdq3b4+hQ4c2+rXUggULoFKpDLfAwECTvg+yfipHOWYMq7uw30e/pKJaqxM5ERER3Q7RBxQbY9++fZg/fz7++9//Ijk5GZs2bcK2bdvwzjvv3HSb2bNno7S01HDLyspqwcRkLSb3D4GvSonc0mp8k3BZ7DhERHQbRCs3np6ekMlkyM/Pr7c8Pz8fPj4+DW4zZ84cTJw4EU8++STCw8Nx//33Y/78+ViwYAH0en2D2ygUCri6uta7Ef2TUi7Di7EdAQBL911EaZVW5ERERNRcopUbe3t7REZGIj4+3rBMr9cjPj4eMTExDW5TWVkJqbR+ZJlMBgA8jZdu2796+6ODlzOuVWrxxf5LYschIqJmEvVrqbi4OCxfvhyrV69GSkoKnnnmGVRUVGDKlCkAgEmTJmH27NmG9UePHo3PP/8c69atQ3p6Onbv3o05c+Zg9OjRhpJD1Fx2Min+c3dnAMDKw+koKOP1k4iIrJGdmC8+fvx4FBYW4s0330ReXh4iIiKwc+dOwyDjzMzMekdq3njjDUgkErzxxhvIyclB27ZtMXr0aLz77rtivQWyMbFdvBAR6IZTWdfw1cF0vDaqi9iRiIjISBKhlX2fo1aroVKpUFpayvE31KC95wow5evjcJDLcPCVYfB0VogdiYio1TPm89uqzpYiaglDO7VFjwAVqrQ6fHUwXew4RERkJKPLzdy5c3H5Mk+VJdslkUjw/B0dAABrEjJQXFEjciIiIjKG0eXmxx9/RPv27TF8+HB899130Gg05shFJKrhXbzQzc8VlTU6rDiUJnYcIiIygtHl5tSpUzh+/Di6deuGmTNnwsfHB8888wyOHz9ujnxEopBIJHh+eN3Rm9VHLuNaJY/eEBFZi2aNuenVqxc+/fRT5ObmYsWKFcjOzsaAAQPQo0cPfPLJJygtLTV1TqIWd1dXb3TxdUW5phYrD3HsDRGRtbitAcWCIECr1aKmpgaCIMDd3R1LlixBYGAg1q9fb6qMRKKoG3tTN+fUqsMZvGoxEZGVaFa5SUpKwrPPPgtfX1+8+OKL6NWrF1JSUrB//35cuHAB7777Lp5//nlTZyVqcSO6+aCTtwvKNLVYdZhHb4iIrIHR5SY8PBz9+vVDeno6VqxYgaysLCxcuBBhYWGGdSZMmIDCwkKTBiUSg1QqwXPD6/7bXnkoHepqHr0hIrJ0Rpebhx56CBkZGdi2bRvGjh3b4LQHnp6eN53IksjajOruiw5ezlBX12L14Qyx4xAR0S0YXW7+HFvzT1VVVXj77bdNEorIkkilEjz759ibIxmoqtGJnIiIiBpjdLmZN28eysvLb1heWVmJefPmmSQUkaW5J9wXgR4OKK6owYakLLHjEBFRI5p15EYikdyw/PTp0/Dw8DBJKCJLYyeTYtqgdgCA5QfTUavj165ERJaqybOCu7u7QyKRQCKRoGPHjvUKjk6nQ3l5OZ5++mmzhCSyBA9GBuLj3eeRWVyJnWfzcG8PP7EjERFRA5pcbhYvXgxBEPDEE09g3rx5UKlUhsfs7e0REhKCmJgYs4QksgQO9jJM7h+Cxb9ewBf703BPuG+DRzGJiEhcTS43kydPBgCEhoaif//+kMvlZgtFZKkmxYRg2f5LOJNTioRLV9E/zFPsSERE9A9NGnOjVqsNP/fq1QtVVVVQq9UN3ohsmYeTPcb3CQQALDvACTWJiCxRk47cuLu748qVK/Dy8oKbm1uDh+L/HGis0/E0WbJtTw5qh2+OXsaB84X4I1eNrn6uYkciIqK/aVK52bNnj+FMqL1795o1EJGlC/RwxD09/PDz6Vx8eeASFj/cS+xIRET0NxJBEASxQ7QktVoNlUqF0tJSuLryX9zUPL/nlOLezw5BJpVg/6yhCHB3FDsSEZFNM+bz2+jr3OzcuROHDh0y3F+6dCkiIiLwyCOPoKSkxPi0RFaou78KA8LaQKcXsOIQJ9QkIrIkRpebWbNmGQYOnzlzBnFxcRg1ahTS09MRFxdn8oBEluqpwe0BAOsSs1BSUSNyGiIi+pPR5SY9PR1du3YFAGzcuBGjR4/G/PnzsXTpUuzYscPkAYks1aAOnuji64oqrQ7/O3ZZ7DhERHSd0eXG3t4elZWVAIBff/0Vd911FwDAw8ODp4JTqyKRSDB9cCgAYE3CZdTUckoGIiJLYHS5GThwIOLi4vDOO+8gMTER99xzDwDg/PnzCAgIMHlAIkt2T7gf2rooUFCmwfYzV8SOQ0REaEa5WbJkCezs7LBhwwZ8/vnn8Pf3BwDs2LEDd999t8kDElkyezspJvULBgCsPJyOVnbyIRGRReKp4ES36Wq5BjEL96CmVo8NT8egT4iH2JGIiGyOMZ/fTZ5b6u/0ej0uXryIgoIC6PX1xxkMHjy4OU9JZLXaOCtwf4Q/1p/IwsrD6Sw3REQiM7rcHD16FI888gguX758wyF4Tr9ArdWUgSFYfyILO3/PQ3ZJJS/qR0QkIqPH3Dz99NPo06cPfv/9dxQXF6OkpMRwKy4uNkdGIovX2ccVA8M8oRfqzpwiIiLxGH3k5sKFC9iwYQPCwsLMkYfIaj0xMASHLhZhbWImZg7vACdFs771JSKi22T0kZvo6GhcvHjRHFmIrNrQjl4I9XRCWXUtNiZnix2HiKjVMvqfls899xxeeukl5OXlITw8HHK5vN7jPXr0MFk4ImsilUowZUAI3vzxLFYdzsBj0cGQSiVixyIianWMPhVcKr3xYI9EIoEgCFYxoJingpM5VWhq0W9BPMqqa7Fich8M7+ItdiQiIptg1lPB09M5AzLRzTgp7DAhKghfHkjDysPpLDdERCIwutwEBwebIweRzZgUE4yvDqbh8MWrSM0rQycfF7EjERG1KkYPKAaAb775BgMGDICfnx8uX6477XXx4sX48ccfTRqOyBoFuDtiRDcfAMCahAxxwxARtUJGl5vPP/8ccXFxGDVqFK5du2YYY+Pm5obFixebOh+RVZoUEwIA2JScg9IqrbhhiIhaGaPLzWeffYbly5fj9ddfh0wmMyzv06cPzpw5Y9JwRNaqXzsPdPJ2QZVWhw1JPC2ciKglGV1u0tPT0atXrxuWKxQKVFRUmCQUkbWTSCSY1L9ufNo3CRnQ61vV/LRERKIyutyEhobi1KlTNyzfuXMnunTpYopMRDZhbIQ/XJR2yLhaiQMXCsWOQ0TUahh9tlRcXBxmzJiB6upqCIKAxMRErF27FgsWLMBXX31ljoxEVslJYYcHIwOx8nA6Vh/JwNBOXmJHIiJqFYwuN08++SQcHBzwxhtvoLKyEo888gj8/PzwySef4OGHHzZHRiKrNTEmGCsPp2Pf+UJkFFUgxNNJ7EhERDavWaeCP/roo7hw4QLKy8uRl5eH7OxsTJ061dTZiKxeqKcThnZqC0EAvj3K2cKJiFpCs8pNUVERTpw4gZSUlHpnTBHRjSZfPy38+xNZqKypFTcMEVErYFS5OXv2LAYPHgxvb29ER0cjKioKXl5euOOOO5CammqujERWbUjHtgjycIS6uhZbTuaKHYeIyOY1udzk5eVhyJAhKCwsxEcffYTt27dj27Zt+OCDD3DlyhUMGjQIBQUF5sxKZJWkUgkmxdSdFr4mIQNGzlVLRERGavKs4K+88gp+/fVXHD58GEqlst5jVVVVGDhwIO666y4sWLDALEFNhbOCkxhKK7XotyAeVVod1k/vh+h2bcSORERkVYz5/G7ykZvdu3fjlVdeuaHYAICDgwNmzZqFXbt2GZ+WqBVQOcoxtpc/AGA155siIjKrJpebtLQ09O7d+6aP9+nTB2lpaSYJRWSLJl+/YvGus/nIV1eLnIaIyHY1udyUlZU1ehjIxcUF5eXlJglFZIs6+7giKsQDOr2AtYmZYschIrJZRp0tVVZWBrVafdMbB0oSNe6x6wOLvzuWCa1OL3IaIiLb1OQrFAuCgI4dOzb6uEQiMUkoIlt1dzcfeDorUFCmwe4/8jEq3FfsSERENqfJ5Wbv3r3mzEHUKtjbSTEhKhCf7bmIbxIus9wQEZlBk8vNkCFDzJmDqNWYEBWEpXsvIiHtKi4WlCHMy0XsSERENqVZ0y8QUfP5uTkgtos3AOCbBM43RURkaiw3RCKYdH2+qY3JOajQcL4pIiJTEr3cLF26FCEhIVAqlYiOjkZiYmKj61+7dg0zZsyAr68vFAoFOnbsiO3bt7dQWiLT6N++Ddp5OqFcU4stp3LEjkNEZFNELTfr169HXFwc5s6di+TkZPTs2RMjRoy46RxVNTU1uPPOO5GRkYENGzYgNTUVy5cvh7+/fwsnJ7o9UqkEj/arOy38m4TLvIwCEZEJGVVutFot7Ozs8Pvvv5vkxT/66CNMmzYNU6ZMQdeuXbFs2TI4Ojpi5cqVDa6/cuVKFBcXY8uWLRgwYABCQkIwZMgQ9OzZ0yR5iFrSA5EBUMqlOJdXhhOXS8SOQ0RkM4wqN3K5HEFBQdDpdLf9wjU1NUhKSkJsbOxfYaRSxMbGIiEhocFtfvrpJ8TExGDGjBnw9vZG9+7dMX/+/EbzaDSaGy42SGQJVA5yjI2oO+rIgcVERKZj9NdSr7/+Ol577TUUFxff1gsXFRVBp9PB29u73nJvb2/k5eU1uE1aWho2bNgAnU6H7du3Y86cOVi0aBH+7//+76avs2DBAqhUKsMtMDDwtnITmdJj17+a2vH7FRSWaUROQ0RkG5p8nZs/LVmyBBcvXoSfnx+Cg4Ph5ORU7/Hk5GSThfsnvV4PLy8vfPnll5DJZIiMjEROTg4++OADzJ07t8FtZs+ejbi4OMN9tVrNgkMWo7u/Cr2C3HAy8xrWH8/Es3d0EDsSEZHVM7rcjB071iQv7OnpCZlMhvz8/HrL8/Pz4ePj0+A2vr6+kMvlkMlkhmVdunRBXl4eampqYG9vf8M2CoUCCoXCJJmJzGFSTDBOZl7Dd8cy8czQMMiknMaEiOh2GF1ubnaExFj29vaIjIxEfHy8oTDp9XrEx8fj2WefbXCbAQMG4LvvvoNer4dUWveN2vnz5+Hr69tgsSGyBiO7++KdrSnILa3GnnMFuLOr9603IiKim2rWqeDXrl3DV199hdmzZxvG3iQnJyMnx7jrdcTFxWH58uVYvXo1UlJS8Mwzz6CiogJTpkwBAEyaNAmzZ882rP/MM8+guLgYM2fOxPnz57Ft2zbMnz8fM2bMaM7bILIISrkMD/Wp+6r0m6McWExEdLuMPnLz22+/ITY2FiqVChkZGZg2bRo8PDywadMmZGZmYs2aNU1+rvHjx6OwsBBvvvkm8vLyEBERgZ07dxoGGWdmZhqO0ABAYGAgdu3ahRdffBE9evSAv78/Zs6ciVdeecXYt0FkUR6NDsIXBy7hwPlCZBRVIMTT6dYbERFRgySCkVcPi42NRe/evfH+++/DxcUFp0+fRrt27XDkyBE88sgjyMjIMFNU01Cr1VCpVCgtLYWrq6vYcYgMpqxKxN7UQkwbFIrX7+kqdhwiIotizOe30V9LHT9+HE899dQNy/39/W96CjcR3drEmLrTwr8/kY1q7e1fS4qIqLUyutwoFIoGL4R3/vx5tG3b1iShiFqjIR29EODugNIqLX4+nSt2HCIiq2V0uRkzZgzefvttaLVaAIBEIkFmZiZeeeUVjBs3zuQBiVoLmVSCR6Prjt58y4HFRETNZnS5WbRoEcrLy+Hl5YWqqioMGTIEYWFhcHFxwbvvvmuOjEStxkN9AmAvk+J0dilOZ10TOw4RkVUy+mwplUqF3bt349ChQ/jtt99QXl6O3r1715sjioiap42zAvf08MXmkzn49uhl9Ax0EzsSEZHVMfpsKWvHs6XI0iVdLsG4z49AYSfFsdeGw82RF6gkIjLm87tJR24+/fRTTJ8+HUqlEp9++mmj6z7//PNNT0pEN+gd5Iauvq7444oaG5Ky8eSgdmJHIiKyKk06chMaGooTJ06gTZs2CA0NvfmTSSRIS0szaUBT45EbsgZrEzMxe9MZhHo6IT5uCKScb4qIWjmTH7lJT09v8GciMo/7Ivwwf1sK0osqcOhiEQZ35GUWiIiayqizpbRaLdq3b4+UlBRz5SEiAI72dhgXGQAAWJPA08KJiIxhVLmRy+Worq42VxYi+pvH+tVd82bPuXxkl1SKnIaIyHoYfZ2bGTNm4L333kNtba058hDRdWFezhgY5gm9APzvWKbYcYiIrIbR17k5fvw44uPj8csvvyA8PBxOTvVnL960aZPJwhG1dhNjgnHoYhHWH8/CzOEdoJTLxI5ERGTxjC43bm5unGaBqIUM7+wFP5USuaXV2H7mCv7VO0DsSEREFs/ocrNq1Spz5CCiBtjJpHi0XzA+2JWK1QmXWW6IiJrA6DE3RNSyxvcNrJtvKusafsu+JnYcIiKLZ/SRm9DQUEgkN7+gmKVfxI/I2ng6KzAq3AdbTuViTcJlfPigm9iRiIgsmtHl5oUXXqh3X6vV4uTJk9i5cydmzZplqlxE9DcTY0Kw5VQufj6di9dHdYG7E+ebIiK6GaPLzcyZMxtcvnTpUpw4ceK2AxHRjXoHuaGbnyvO5qrx/YksPDWkvdiRiIgslsnG3IwcORIbN2401dMR0d9IJBJMiqm7qN+3xy5Dp7/llHBERK2WycrNhg0b4OHhYaqnI6J/GNPTHyoHObKKq7D/fIHYcYiILJbRX0v16tWr3oBiQRCQl5eHwsJC/Pe//zVpOCL6i4O9DA9GBuCrQ+lYfeQy7ujsLXYkIiKLZHS5GTt2bL37UqkUbdu2xdChQ9G5c2dT5SKiBkyMCcaKw+nYf74Q6UUVCPV0uvVGREStjNHlZu7cuebIQURNENzGCcM6eWHPuQKsScjA3NHdxI5ERGRxjC43AKDT6bBlyxakpKQAALp164YxY8ZAJuO8N0TmNrl/CPacK8APJ7Lx0l2d4Kxo1q8xEZHNMnpA8cWLF9GlSxdMmjQJmzZtwqZNm/DYY4+hW7duuHTpkjkyEtHfDArzRLu2TijX1GJTcrbYcYiILI7R5eb5559H+/btkZWVheTkZCQnJyMzMxOhoaF4/vnnzZGRiP5GKpVgckwIAODrIxnQ87RwIqJ6jC43+/fvx/vvv1/vtO82bdpg4cKF2L9/v0nDEVHDxkUGwFlhh7TCChy8WCR2HCIii2J0uVEoFCgrK7theXl5OezteUl4opbgrLDDg33qZghffSRD3DBERBbG6HJz7733Yvr06Th27BgEQYAgCDh69CiefvppjBkzxhwZiagBk65/NbU3tQAZRRXihiEisiBGl5tPP/0U7du3R0xMDJRKJZRKJQYMGICwsDB88skn5shIRA0I9XTCsE5tIQjAmoTLYschIrIYRp9D6ubmhh9//BEXLlzAuXPnAABdunRBWFiYycMRUeMm9w/B3tRC/HAiCy/d1RFOPC2ciKjpR27efPNNVFZWGu57enpi9OjRGD16NIsNkUgGd2iLdp5OKONp4UREBk0uN++++y7Ky8sN94ODg5GWlmaWUETUNFLpX7OF87RwIqI6TS43giA0ep+IxPHnaeGXCitwiKeFExEZP6CYiCyLi1KOByLrTgtfdThd5DREROJrcrmRSCQoKyuDWq1GaWkpJBIJysvLoVar692IqOU93j8EEgmwN7UQFwvKb70BEZENM+prqY4dO8Ld3R0eHh4oLy9Hr1694O7uDnd3d7i5ucHd3d2cWYnoJkI8nTC8szcAHr0hImryeaN79+41Zw4iuk1TB4bi15R8bEzOxst3dYK7E68YTkStU5PLzZAhQ8yZg4huU792Hujm54qzuWp8l5iJGcN4iQYiap04oJjIRkgkEkwdGAqgbr6pmlq9yImIiMTBckNkQ+7t4QcvFwUKyjTY+luu2HGIiETBckNkQ+ztpJjcPwQAsOJQOq9HRUStEssNkY15JCoISrkUZ3PVOJpWLHYcIqIWx3JDZGPcnewxrnfdRf1WHOJp4UTU+jTpbKl//etfTX7CTZs2NTsMEZnGEwND8b9jmYg/l4/0ogqEejqJHYmIqMU06ciNSqUy3FxdXREfH48TJ04YHk9KSkJ8fDxUKpXZghJR07Vv64w7OntBEHhRPyJqfZp05GbVqlWGn1955RU89NBDWLZsGWQyGQBAp9Ph3//+N1xdXc2TkoiMNnVgKPacK8APJ7Lx0p2doHKUix2JiKhFGD3mZuXKlXj55ZcNxQYAZDIZ4uLisHLlSpOGI6Lm69++DTr7uKBKq8O3xy6LHYeIqMUYXW5qa2tx7ty5G5afO3cOej0vGkZkKSQSCaYPbgcAWHU4A9VanciJiIhaRpOnX/jTlClTMHXqVFy6dAlRUVEAgGPHjmHhwoWYMmWKyQMSUfON7umHD3elIre0GpuSc/BIdJDYkYiIzM7ocvPhhx/Cx8cHixYtwpUrVwAAvr6+mDVrFl566SWTBySi5pPLpJg6qB3e2foHlh9Mw/i+gZBJJWLHIiIyK4lwG5cwVavVAGBVA4nVajVUKhVKS0utKjdRc1VoatF/4R6UVmnx+aO9MTLcV+xIRERGM+bzu1kX8autrcWvv/6KtWvXQiKp+1dgbm4uysvLm/N0RGRGTgo7TIoJBgAs23+JUzIQkc0zutxcvnwZ4eHhuO+++zBjxgwUFhYCAN577z28/PLLJg9IRLdvcv8QKOykOJ1dyikZiMjmGV1uZs6ciT59+qCkpAQODg6G5ffffz/i4+NNGo6ITMPTWYEH+9RNybBs/yWR0xARmZfR5ebgwYN44403YG9vX295SEgIcnJyTBaMiExr2qB2kEqA/ecLkXJFLXYcIiKzMbrc6PV66HQ3Xi8jOzsbLi4uzQqxdOlShISEQKlUIjo6GomJiU3abt26dZBIJBg7dmyzXpeoNQlu42QYTPwFj94QkQ0zutzcddddWLx4seG+RCJBeXk55s6di1GjRhkdYP369YiLi8PcuXORnJyMnj17YsSIESgoKGh0u4yMDLz88ssYNGiQ0a9J1Fo9Pbg9AODn364gu6RS5DREROZhdLlZtGgRDh8+jK5du6K6uhqPPPKI4Sup9957z+gAH330EaZNm4YpU6aga9euWLZsGRwdHRudykGn0+HRRx/FvHnz0K5dO6Nfk6i1Cg9QYUBYG+j0Ar46yAk1icg2GV1uAgICcPr0abz++ut48cUX0atXLyxcuBAnT56El5eXUc9VU1ODpKQkxMbG/hVIKkVsbCwSEhJuut3bb78NLy8vTJ069ZavodFooFar692IWrOnh9QdvVl/PAvFFTUipyEiMj2jr1AMAHZ2dnj00Ufx6KOP3taLFxUVQafTwdvbu95yb2/vBuevAoBDhw5hxYoVOHXqVJNeY8GCBZg3b95t5SSyJQPDPBHur8KZnFKsOJSGWSM6ix2JiMikjD5yI5PJMGzYMBQX179WRn5+fr2Zws2hrKwMEydOxPLly+Hp6dmkbWbPno3S0lLDLSsry6wZiSydRCLBs3eEAQBWH7mM0kqtyImIiEzL6CM3giBAo9GgT58++Pnnn9GtW7d6jxnD09MTMpkM+fn59Zbn5+fDx8fnhvUvXbqEjIwMjB492rDsz5nI7ezskJqaivbt29fbRqFQQKFQGJWLyNbd2cUbnX1ccC6vDKuOpOOF2I5iRyIiMhmjj9xIJBJs3LgRo0ePRkxMDH788cd6jxnD3t4ekZGR9S7+p9frER8fj5iYmBvW79y5M86cOYNTp04ZbmPGjMGwYcNw6tQpBAYGGvt2iFolqfSvozcrD6WjrJpHb4jIdjTryI1MJsMnn3yCbt26Yfz48XjjjTfw5JNPNitAXFwcJk+ejD59+iAqKgqLFy9GRUUFpkyZAgCYNGkS/P39sWDBAiiVSnTv3r3e9m5ubgBww3IiatzI7r5o3/Y8LhVWYE3CZcwYFiZ2JCIik2jWgOI/TZ8+HR06dMCDDz6IAwcONOs5xo8fj8LCQrz55pvIy8tDREQEdu7caRhknJmZCam0WfN7ElEjZNeP3ry4/jS+OpiGx/uHwElxW/9LICKyCBLByIEyoaGhOHHiBNq0aWNYdvHiRYwePRrnz59v8OrFlsSYKdOJbF2tTo/hH+3H5auVeG1UZ0wf3P7WGxERicCYz2+jD4mkp6fXKzYAEBYWhpMnTyItLc3YpyMiEdnJpJgxtO7rqC8PpKNaa9n/OCEiagqTfd+jVCoRHBxsqqcjohZyf29/+Ls5oKhcg7WJmWLHISK6bU0qNx4eHigqKgIAuLu7w8PD46Y3IrIucpkUzwyt+zrqi/1p0NTy6A0RWbcmjR78+OOPDTN+/33STCKyDQ/2CcCSPReRp67GDyey8Vg/HoUlIutl9IBia8cBxUQN+/pwOt76+Q/4qZTYO2soFHbmveI4EZExjPn8btKRG2Mmm2RhILJOD0cFYdn+NOSWVmPtsUw8PiBU7EhERM3SpHLj5uZ2y6sPC4IAiURi8aeCE1HDlHIZnr0jDG9s+R1L9l7CQ30D4WjP694QkfVp0v+59u7da+4cRGQBHuoTiC8OXEJWcRVWH7lsGGhMRGRNOOaGiOrZmJSNl344DZWDHAdfGQZXpVzsSEREph9z05DKykpkZmaipqam3vIePXo09ymJyAKM7eWPz/dfwsWCcnx1MB1xd3LGcCKyLkaXm8LCQkyZMgU7duxo8HGOuSGybjKpBHF3dsS//5eMFdfnnPJwshc7FhFRkxl9heIXXngB165dw7Fjx+Dg4ICdO3di9erV6NChA3766SdzZCSiFnZ3Nx9083NFRY0Oy/ZfEjsOEZFRjC43e/bswUcffYQ+ffpAKpUiODgYjz32GN5//30sWLDAHBmJqIVJpRK8fFcnAMDqIxnIV1eLnIiIqOmMLjcVFRXw8vICUDcVQ2FhIQAgPDwcycnJpk1HRKIZ2qktIoPdoanVY8mei2LHISJqMqPLTadOnZCamgoA6NmzJ7744gvk5ORg2bJl8PX1NXlAIhKHRPLX0Zt1xzORVVwpciIioqYxutzMnDkTV65cAQDMnTsXO3bsQFBQED799FPMnz/f5AGJSDwx7dtgYJgntDoBH+0+L3YcIqImue3r3FRWVuLcuXMICgqCp6enqXKZDa9zQ2ScM9mlGL3kEADg52cHIjxAJXIiImqNjPn8NvrIzT85Ojqid+/eVlFsiMh44QEq3N/LHwDw7vY/0Mqu+0lEVsjo69wIgoANGzZg7969KCgogF6vr/f4pk2bTBaOiCzDyyM6YduZKziaVoxfUwpwZ1dvsSMREd1Us65zM3HiRKSnp8PZ2RkqlarejYhsj7+bA6YOrJslfMGOFGh1+ltsQUQkHqOP3HzzzTfYtGkTRo0aZY48RGSh/j20Pb4/noW0wgqsS8zExJgQsSMRETXI6CM3KpUK7dq1M0cWIrJgLko5XojtAAD4+NcLUFdrRU5ERNQwo8vNW2+9hXnz5qGqqsoceYjIgj0cFYR2bZ1QXFGDz/dxWgYiskxGl5uHHnoIJSUl8PLyQnh4OHr37l3vRkS2Sy6T4rWRXQAAKw6lI+ca/5FDRJbH6DE3kydPRlJSEh577DF4e3tDIpGYIxcRWajhXbzQr50HjqYV48Ndqfh4fITYkYiI6jH6In5OTk7YtWsXBg4caK5MZsWL+BHdvr9f2G/zv/ujV5C7yImIyNaZ9SJ+gYGBLAVErVx4gArjegcAAOb+dBY6PS/sR0SWw+hys2jRIvznP/9BRkaGGeIQkbV4dWRnuCjs8Ft2KdYdzxQ7DhGRgdFfS7m7u6OyshK1tbVwdHSEXC6v93hxcbFJA5oav5YiMp1Vh9Mx7+c/4OYox96XhsLdyV7sSERko4z5/DZ6QPHixYubm4uIbMzEfsFYfzwL5/LK8P6uVCz4V7jYkYiIjCs3Wq0W+/fvx5w5cxAaGmquTERkJexkUrx9X3c89EUC1h3PxMN9A9Ez0E3sWETUyhk15kYul2Pjxo3mykJEVigq1AP/6uUPQQDe/PF36Dm4mIhEZvSA4rFjx2LLli1miEJE1urVUXWDi09nl2L9iSyx4xBRK2f0mJsOHTrg7bffxuHDhxEZGQknJ6d6jz///PMmC0dE1sHLRYkX7uyId7b+gfd2nsPd3Xw4uJiIRGP02VKNjbWRSCRIS0u77VDmxLOliMyjVqfHPZ8eQmp+GSZEBXFwMRGZlFnPlkpPT292MCKyXXWDi7th/JdHsTYxE/f38kdUqIfYsYioFTJ6zM3fCYIAIw/8EJENi27XBuP7BAIAXt34G6q1OpETEVFr1Kxys2bNGoSHh8PBwQEODg7o0aMHvvnmG1NnIyIr9No9XeDlokBaUQU+jb8gdhwiaoWMLjcfffQRnnnmGYwaNQrff/89vv/+e9x99914+umn8fHHH5sjIxFZEZWDHO+M7Q4A+OJAGs7mloqciIham2YNKJ43bx4mTZpUb/nq1avx1ltvWfyYHA4oJmoZM/6XjG1nrqCbnyt+nDEAdrLb+haciFo5s84KfuXKFfTv3/+G5f3798eVK1eMfToislFvjekGlYMcZ3PVWH7Qsv/RQ0S2xehyExYWhu+///6G5evXr0eHDh1MEoqIrF9bFwXevLcrAODjX88jrbBc5ERE1FoYfSr4vHnzMH78eBw4cAADBgwAABw+fBjx8fENlh4iar3+1dsfP57OxYHzhXh14xmsm94PUqlE7FhEZOOMPnIzbtw4HDt2DJ6entiyZQu2bNkCT09PJCYm4v777zdHRiKyUhKJBPPv7w5HexkSM4rx7bHLYkciolbA6AHF1o4Diola3uojGZj701ko5VJse34Q2rd1FjsSEVkZsw4oJiIy1sR+wRgY5olqrR4vrDuFmlq92JGIyIY1udxIpVLIZLJGb3Z2Rg/hIaJWQCqVYNFDPeHmKMeZnFJ8En9e7EhEZMOa3EY2b95808cSEhLw6aefQq/nv8aIqGHerkrMvz8c//5fMv677xKGdPTi3FNEZBa3NeYmNTUVr776Kn7++Wc8+uijePvttxEcHGzKfCbHMTdE4nr5h9PYkJQNfzcH7HhhEFyVcrEjEZEVMPuYm9zcXEybNg3h4eGora3FqVOnsHr1aosvNkQkvrfGdEOQhyNyrlXhrR/Pih2HiGyQUeWmtLQUr7zyCsLCwnD27FnEx8fj559/Rvfu3c2Vj4hsjLPCDh+P7wmpBNh0Mgc/n84VOxIR2Zgml5v3338f7dq1w9atW7F27VocOXIEgwYNMmc2IrJRkcEeeHZYGADg9c1nkHOtSuRERGRLmjzmRiqVwsHBAbGxsZDJZDddb9OmTSYLZw4cc0NkGbQ6PR5YloDTWdcQEeiG75+Kgb0dr05BRA0z5vO7yWdLTZo0CRIJL5tORKYhl0mxZEIv3PvZIZzKuob521Pw1phuYsciIhvAKxQTkajiU/IxdfUJAMCnE3phTE8/kRMRkSXiFYqJyGoM7+KNfw9tDwB4deNvuFhQJnIiIrJ2FlFuli5dipCQECiVSkRHRyMxMfGm6y5fvhyDBg2Cu7s73N3dERsb2+j6RGT54u7siJh2bVBZo8PT3yajQlMrdiQismKil5v169cjLi4Oc+fORXJyMnr27IkRI0agoKCgwfX37duHCRMmYO/evUhISEBgYCDuuusu5OTktHByIjIVO5kUn07oBS8XBS4WlGP2pjNoZd+YE5EJiT7mJjo6Gn379sWSJUsAAHq9HoGBgXjuuefw6quv3nJ7nU4Hd3d3LFmyBJMmTbrl+hxzQ2S5jmcU4+Evj0KnF/DOfd0wMSZE7EhEZCGsZsxNTU0NkpKSEBsba1gmlUoRGxuLhISEJj1HZWUltFotPDwanqNGo9FArVbXuxGRZeob4oFX7+4MAHh76x84mnZV5EREZI1ELTdFRUXQ6XTw9vaut9zb2xt5eXlNeo5XXnkFfn5+9QrS3y1YsAAqlcpwCwwMvO3cRGQ+Tw4KxT09fKHVCXj62yRkFFWIHYmIrIzoY25ux8KFC7Fu3Tps3rwZSqWywXVmz56N0tJSwy0rK6uFUxKRMSQSCT58oCd6BqhwrVKLqauPo7RKK3YsIrIiopYbT09PyGQy5Ofn11uen58PHx+fRrf98MMPsXDhQvzyyy/o0aPHTddTKBRwdXWtdyMiy+ZgL8PySX3gq1LiUmEFnv0uGVqdXuxYRGQlRC039vb2iIyMRHx8vGGZXq9HfHw8YmJibrrd+++/j3feeQc7d+5Enz59WiIqEbUwL1clvprcB472Mhy8UIR5P5/lGVRE1CSify0VFxeH5cuXY/Xq1UhJScEzzzyDiooKTJkyBUDdtA+zZ882rP/ee+9hzpw5WLlyJUJCQpCXl4e8vDyUl5eL9RaIyEy6+amweHwEJBLg26OZWH0kQ+xIRGQFRC8348ePx4cffog333wTEREROHXqFHbu3GkYZJyZmYkrV64Y1v/8889RU1ODBx54AL6+vobbhx9+KNZbICIzuqubT70zqPamNnwNLCKiP4l+nZuWxuvcEFkfQRDwysbf8P2JbDjIZfjftGj0DnIXOxYRtSCruc4NEVFTSCQS/N/YcAzq4IkqrQ5TVh1Hah7noCKihrHcEJFVsLeT4ouJkegV5IbSKi0mrjiGrOJKsWMRkQViuSEiq+Fob4dVj/dFJ28XFJRpMHHFMRSWacSORUQWhuWGiKyKm6M91kyNQqCHAzKuVmLSykRe5I+I6mG5ISKr4+2qxDdPRMPTWYGUK2o8ufo4qmp0YsciIgvBckNEVinE0wlrnoiCi9IOxzNKMHX1cVTW1Iodi4gsAMsNEVmtrn6u+HpKXzjZy3Dk0lU8vuo4KjQsOEStHcsNEVm1yGAPrJkaDReFHRLTizF5ZSLKqjkGh6g1Y7khIqsXGeyOb56MhqvSDicul2DSykSoWXCIWi2WGyKyCRGBbvhuWj+oHOQ4mXkNE786htJKFhyi1ojlhohsRnd/FdZO6wd3RzlOZ5fika+Ooqic18Eham1YbojIpnT1c8Xa6f3QxskeZ3PVGPf5EWQUVYgdi4haEMsNEdmczj6u+OHpGAR6OODy1UqM+/wITmVdEzsWEbUQlhsisknt2jpj4zP90d3fFVcrajDhy6PYcy5f7FhE1AJYbojIZnm5KLFuegwGd2yLKq0O09YkYf3xTLFjEZGZsdwQkU1zVthhxeQ+GNc7ADq9gFc2nsFHu89DrxfEjkZEZsJyQ0Q2Ty6T4sMHe+DZYWEAgE/jL+Dpb5NQzqsZE9kklhsiahUkEgleHtEJ7z/QA/YyKX75Ix/3Lz2MdJ5JRWRzWG6IqFV5qE8g1j/VD96uClwoKMeYJYewN7VA7FhEZEIsN0TU6vQKcsfPzw1EZLA7yqpr8cTXx/HffRchCByHQ2QLWG6IqFXyclFi7bR+mBAVBEEA3t+ZimlrklBSUSN2NCK6TSw3RNRq2dtJseBf4Xj3/u6wl0nxa0o+Rn5yEAmXroodjYhuA8sNEbV6j0YHY9O/+6NdWyfkqavxyFdHseiXVNTq9GJHI6JmYLkhIkLdpJtbnxuIh/oEQBCAz/ZcxPgvjyKruFLsaERkJJYbIqLrHO3t8P4DPfHZhF5wUdgh6XIJRn16EN+fyOJgYyIrwnJDRPQPo3v6YfvMQegd5Iay6lr8Z8NvmLQyEdklPIpDZA1YboiIGhDo4Yjvn4rB7JGdobCT4uCFIoz4+ADWJGRw6gYiC8dyQ0R0E3YyKZ4a0h47Zg5C3xB3VNTo8OaPZ/Hwl0d5ZWMiC8ZyQ0R0C+3aOmP99BjMG9MNjvYyJGYUY8THB/DBrnOorOH8VESWhuWGiKgJpFIJJvcPwa4XBmNwx7ao0emxdO8lDF+0Hz+fzuWAYyILwnJDRGSEQA9HrJ7SF19MjESAuwOulFbjubUnMWH5UZzLU4sdj4gASIRW9s8NtVoNlUqF0tJSuLq6ih2HiKxYtVaHL/an4b/7LkJTq4dMKsHDfQPx/PAO8HZVih2PyKYY8/nNckNEdJuySyrx7rYU7Pg9DwCglEvxeP9QPDOkPVSOcpHTEdkGlptGsNwQkbkcTbuK93eeQ3LmNQCAq9IOTw9tjyn9Q+FgLxM3HJGVY7lpBMsNEZmTIAj4NaUAH+5KRWp+GQCgrYsCTw1uhwlRQXBS2ImckMg6sdw0guWGiFqCTi/gp9M5WPTLeWSXVAEA3BzlmNI/FJP7B8PN0V7khETWheWmESw3RNSSamr12HwyG5/vu4SMq3XTNzjZy/BYv2BMHRgKLw48JmoSlptGsNwQkRh0egHbz1zB0r0XcS6v7usqe5kU9/TwxeP9Q9Az0E3cgEQWjuWmESw3RCQmQRCwN7UAS/deQtLlEsPyiEA3PN4/BKPCfWFvx0uQEf0Ty00jWG6IyFKczrqG1UcysPW3K6jR6QEAns4KPNw3EA/2CUBwGyeRExJZDpabRrDcEJGlKSzTYG1iJr49ehkFZRrD8n7tPPBQn0CM7O7LU8mp1WO5aQTLDRFZqppaPX75Iw/fn8jGwQuF+PP/zs4KO4zu6YuxEf7oG+IBqVQiblAiEbDcNILlhoisQc61KmxKysb3SVnIKq4yLPdxVeKeHr4Y3dMPPQNUkEhYdKh1YLlpBMsNEVkTvV7A0fSr2JScg12/56FMU2t4LNDDAfeE++Gubt6ICHDjER2yaSw3jWC5ISJrpanV4cD5Ivx8Ohe7/8hHlVZneKytiwKxXbxwZ1dv9G/vCaWcY3TItrDcNILlhohsQWVNLfacK8DO3/OwP7Ww3hEdR3sZBoR5YnAHTwzu2JZnXZFNYLlpBMsNEdmamlo9jqZdxe4/8vFrSj6ulFbXezzIwxGDO3piUIe26NeuDVQOnKmcrA/LTSNYbojIlgmCgLO5auw/X4gD5wuRdLkEtfq//jcvkQBdfFwR3c4D0aFtEB3qAXcnznNFlo/lphEsN0TUmpRranH00lUcvFCIgxeKkFZUccM6Hb2dERnsjl6B7ugV5Ib2bZ05OJksDstNI1huiKg1K1BX41h6MY6lX8WxtGJcKCi/YR0XpR0iAt0QEeiG7v4qhPur4KtS8rRzEhXLTSNYboiI/lJUrsGJjBKczCrBycxr+C37Gqq1+hvWa+Nkj27+KnT3c0UXX1d09nFBqKcT7GScB4taBstNI1huiIhuTqvTIzWvDMmZJfgtuxS/55TiQkE5dPobPyrsZVK093JGZx8XdPJxQQcvZ7Rv64wAdweWHjI5Yz6/7VooExERWQG5TIru/ip091cZllVrdUjNK8OZnFKczS1FypUynM8vQ2WNDilX1Ei5oq73HPYyKUI8HRHm5Yx2ns4IbuOIEE8nBLdxRFtnBb/eIrNjuSEiokYp5TL0DHRDz0A3wzK9XkB2SRXO5amRmleG1PwyXCqsQFphOTS1epzPL8f5/BvH8zjayxDk4YggD0cEuDsi0MMBAe6OCHB3QIC7A1yUPE2dbh+/liIiIpPR6wXkXKvCxcJyXCooR1pRBTKvViLjagVyr1WhgW+36nFR2MHXTQkflQP8VEr4qhzgq1LCy1UBLxclvF0VcHe059lcrRDH3DSC5YaISBw1tXpkl9QVneySKmSXVCGruPL6z5UoqdQ26XnkMgnaOivQ1kUBT+frNxd7w89tnOzh7mSPNk72cHO0h70dx//YAo65ISIii2NvJ0W7ts5o19a5wccrNLW4UlqFK6XVuHKtGrmlVcgrrUZuaTUK1NUoLNPgakUNtDoBudeXN4WLwg7uTvZwd5RD5WgPNwc53BzlcHOQw9VBDtX1P12Vf/5sBxelHM4KO8h4hMgqWUS5Wbp0KT744APk5eWhZ8+e+OyzzxAVFXXT9X/44QfMmTMHGRkZ6NChA9577z2MGjWqBRMTEZGpOSnsEOblgjAvl5uuU1OrR1G5Bvl/KztFZRoUlWtQVF6DwnINSipqUFJZg+KKGugFoExTizJNLTKLm5HJXgZn5V9lx1lhByeFDE6Gn+v+dJDL4Ggvg6PCDo5yGRwVMjja1y13kMvgYH/9JpexMLUA0cvN+vXrERcXh2XLliE6OhqLFy/GiBEjkJqaCi8vrxvWP3LkCCZMmIAFCxbg3nvvxXfffYexY8ciOTkZ3bt3F+EdEBFRS7G3k8LPzQF+bg63XFevF6Cu1uJqRQ1KKmpwrVKLa1VaXKusQWmV1nC/rFqL0iot1FVaqKtrUVqlRU1t3bV+Kmp0qKjRIV+tMdl7kMskUMpl129SKO3qflbYSaGQS6Gwq1uusJPBXiaFvZ0UCrs//5TB/vrP9jIJ7O2kkF9fRy6Twl5W96dcJoHd9ft2Mkm9ZXJp3Z92MgnspBLYSeses6Wz2EQfcxMdHY2+fftiyZIlAAC9Xo/AwEA899xzePXVV29Yf/z48aioqMDWrVsNy/r164eIiAgsW7bslq/HMTdERHQrmlodyqtrUVZdi3JNLdTVWpRX16KiphblGh0qNLWo1PztZ60OVTW1qNDoUFlTi4oaHapqdKjW6lB1/WbpI1ylEtSVHqkEMqnk+p9/uy+TQCap+1n2t3Wk0rrlUsM2EnTwcsGbo7uaNJ/VjLmpqalBUlISZs+ebVgmlUoRGxuLhISEBrdJSEhAXFxcvWUjRozAli1bGlxfo9FAo/mrcavV6gbXIyIi+pPCTgaFswxtnBUmeT5BEKCp1dcVnlodqrV//7nuVlOrh6ZWD41WD02tDppa/V/LdXrD4zXXb9rry2r+9metToBWp79+q/9zrU4Prb7uz4bOWtMLdV/71Zjg/VZoak3wLM0narkpKiqCTqeDt7d3veXe3t44d+5cg9vk5eU1uH5eXl6D6y9YsADz5s0zTWAiIqJmkEj++irKEuj0Amr1euj0ArQ6oe7+9fKj19eVorp1hOvr6KEXBNTqBOgE4a/Hrt/X6/9artMLos80L/qYG3ObPXt2vSM9arUagYGBIiYiIiISV93XSpZRtMxB1HLj6ekJmUyG/Pz8esvz8/Ph4+PT4DY+Pj5Gra9QKKBQmOawIhEREVk+Ua9sZG9vj8jISMTHxxuW6fV6xMfHIyYmpsFtYmJi6q0PALt3777p+kRERNS6iP61VFxcHCZPnow+ffogKioKixcvRkVFBaZMmQIAmDRpEvz9/bFgwQIAwMyZMzFkyBAsWrQI99xzD9atW4cTJ07gyy+/FPNtEBERkYUQvdyMHz8ehYWFePPNN5GXl4eIiAjs3LnTMGg4MzMTUulfB5j69++P7777Dm+88QZee+01dOjQAVu2bOE1boiIiAiABVznpqXxOjdERETWx5jPb84mRkRERDaF5YaIiIhsCssNERER2RSWGyIiIrIpLDdERERkU1huiIiIyKaw3BAREZFNYbkhIiIimyL6FYpb2p/XLFSr1SInISIioqb683O7KdcebnXlpqysDAAQGBgochIiIiIyVllZGVQqVaPrtLrpF/R6PXJzc+Hi4gKJRGJY3rdvXxw/frzBbW722D+Xq9VqBAYGIisrS/SpHRp7Py39fMZs25R1m7OvGnuM+9H023I/1mfq/Xg7z9mS+7Gxx61xPwKW8zvZGvejIAgoKyuDn59fvTknG9LqjtxIpVIEBATcsFwmk930L/xmj91suaurq+i/hI29n5Z+PmO2bcq6zdlXjT3G/Wj6bbkf6zP1fryd52zJ/djY49a4HwHL+Z1srfvxVkds/sQBxdfNmDHD6Mca20Zsps52O89nzLZNWbc5+6qxx7gfTb8t92N95sjW3Odsyf3Y2OPWuB8By/md5H5sXKv7WsqcOOO4beB+tA3cj7aB+9E2tPR+5JEbE1IoFJg7dy4UCoXYUeg2cD/aBu5H28D9aBtaej/yyA0RERHZFB65ISIiIpvCckNEREQ2heWGiIiIbArLDREREdkUlhsiIiKyKSw3IsjKysLQoUPRtWtX9OjRAz/88IPYkaiZ7r//fri7u+OBBx4QOwoZYevWrejUqRM6dOiAr776Suw4dBv4O2j9zPGZyFPBRXDlyhXk5+cjIiICeXl5iIyMxPnz5+Hk5CR2NDLSvn37UFZWhtWrV2PDhg1ix6EmqK2tRdeuXbF3716oVCpERkbiyJEjaNOmjdjRqBn4O2j9zPGZyCM3IvD19UVERAQAwMfHB56eniguLhY3FDXL0KFD4eLiInYMMkJiYiK6desGf39/ODs7Y+TIkfjll1/EjkXNxN9B62eOz0SWmwYcOHAAo0ePhp+fHyQSCbZs2XLDOkuXLkVISAiUSiWio6ORmJjYrNdKSkqCTqdDYGDgbaamf2rJ/Ugt53b3a25uLvz9/Q33/f39kZOT0xLR6R/4O2obTLkfTfWZyHLTgIqKCvTs2RNLly5t8PH169cjLi4Oc+fORXJyMnr27IkRI0agoKDAsE5ERAS6d+9+wy03N9ewTnFxMSZNmoQvv/zS7O+pNWqp/UgtyxT7lSwD96VtMNV+NOlnokCNAiBs3ry53rKoqChhxowZhvs6nU7w8/MTFixY0OTnra6uFgYNGiSsWbPGVFGpEebaj4IgCHv37hXGjRtniphkpObs18OHDwtjx441PD5z5kzhf//7X4vkpZu7nd9R/g5ajubuR1N/JvLIjZFqamqQlJSE2NhYwzKpVIrY2FgkJCQ06TkEQcDjjz+OO+64AxMnTjRXVGqEKfYjWZ6m7NeoqCj8/vvvyMnJQXl5OXbs2IERI0aIFZlugr+jtqEp+9Ecn4ksN0YqKiqCTqeDt7d3veXe3t7Iy8tr0nMcPnwY69evx5YtWxAREYGIiAicOXPGHHHpJkyxHwEgNjYWDz74ILZv346AgAD+T1dkTdmvdnZ2WLRoEYYNG4aIiAi89NJLPFPKAjX1d5S/g5atKfvRHJ+Jdre1NTXLwIEDodfrxY5BJvDrr7+KHYGaYcyYMRgzZozYMcgE+Dto/czxmcgjN0by9PSETCZDfn5+veX5+fnw8fERKRUZi/vRNnG/2g7uS9sg1n5kuTGSvb09IiMjER8fb1im1+sRHx+PmJgYEZORMbgfbRP3q+3gvrQNYu1Hfi3VgPLycly8eNFwPz09HadOnYKHhweCgoIQFxeHyZMno0+fPoiKisLixYtRUVGBKVOmiJia/on70TZxv9oO7kvbYJH70STnXNmYvXv3CgBuuE2ePNmwzmeffSYEBQUJ9vb2QlRUlHD06FHxAlODuB9tE/er7eC+tA2WuB85txQRERHZFI65ISIiIpvCckNEREQ2heWGiIiIbArLDREREdkUlhsiIiKyKSw3REREZFNYboiIiMimsNwQERGRTWG5ISLRDB06FC+88MJtPYdEIoFEIoGbm5th2VtvvYWIiIjbet6W8NZbbxnyL168WOw4RDaD5YaIGlRYWIhnnnkGQUFBUCgU8PHxwYgRI3D48GGxo91g1apVOH/+fLO2nTJlCt544w0AfxWlo0eP1ltHo9GgTZs2kEgk2Ldv3+3GNXj55Zdx5coVBAQEmOw5iYgTZxLRTYwbNw41NTVYvXo12rVrh/z8fMTHx+Pq1atiR7uBm5sbvLy8jN5Op9Nh69at2LZtm2FZYGAgVq1ahX79+hmWbd68Gc7OziguLjZJ3j85OzvD2dkZMpnMpM9L1NrxyA0R3eDatWs4ePAg3nvvPQwbNgzBwcGIiorC7NmzMWbMmHrrPfXUU/D29oZSqUT37t2xdetWAMDVq1cxYcIE+Pv7w9HREeHh4Vi7dm2jr6vRaPDyyy/D398fTk5OiI6Ovq0jJd988w1CQkKgUqnw8MMPo6ysrN7jR44cgVwuR9++fQ3LJk+ejHXr1qGqqsqwbOXKlZg8eXK9bTMyMiCRSLBu3Tr079/f8P73799fb72zZ8/i3nvvhaurK1xcXDBo0CBcunSp2e+JiG6N5YaIbvDnEYUtW7ZAo9E0uI5er8fIkSNx+PBhfPvtt/jjjz+wcOFCw1GI6upqREZGYtu2bfj9998xffp0TJw4EYmJiTd93WeffRYJCQlYt24dfvvtNzz44IO4++67ceHCBaPfw6VLl7BlyxZs3boVW7duxf79+7Fw4cJ66/z0008YPXo0JBKJYVlkZCRCQkKwceNGAEBmZiYOHDiAiRMnNvg6s2bNwksvvYSTJ08iJiYGo0ePNhzdysnJweDBg6FQKLBnzx4kJSXhiSeeQG1trdHvh4iMYNY5x4nIam3YsEFwd3cXlEql0L9/f2H27NnC6dOnDY/v2rVLkEqlQmpqapOf85577hFeeuklw/0hQ4YIM2fOFARBEC5fvizIZDIhJyen3jbDhw8XZs+efdPnBCBs3ry53rK5c+cKjo6OglqtNiybNWuWEB0dXW+9Dh06CFu3br3huRYvXiwMGzZMEARBmDdvnnD//fcLJSUlAgBh7969giAIQnp6ugBAWLhwoWF7rVYrBAQECO+9954gCIIwe/ZsITQ0VKipqWn07yU4OFj4+OOPG12HiJqOR26IqEHjxo1Dbm4ufvrpJ9x9993Yt28fevfuja+//hoAcOrUKQQEBKBjx44Nbq/T6fDOO+8gPDwcHh4ecHZ2xq5du5CZmdng+mfOnIFOp0PHjh0NR46cnZ2xf//+Zn2NExISAhcXF8N9X19fFBQUGO6npKQgNzcXw4cPv2Hbxx57DAkJCUhLS8PXX3+NJ5544qavExMTY/jZzs4Offr0QUpKCoC6v6NBgwZBLpcbnZ+Imo8DionoppRKJe68807ceeedmDNnDp588knMnTsXjz/+OBwcHBrd9oMPPsAnn3yCxYsXIzw8HE5OTnjhhRdQU1PT4Prl5eWQyWRISkq6YYCts7Oz0dn/WSgkEgn0er3h/k8//YQ777wTSqXyhm3btGmDe++9F1OnTkV1dTVGjhx5w3idprjV3xERmQeP3BBRk3Xt2hUVFRUAgB49eiA7O/ump2AfPnwY9913Hx577DH07NkT7dq1a/R07V69ekGn06GgoABhYWH1bj4+PiZ/Lz/++CPuu+++mz7+xBNPYN++fZg0aVKjZzP9/bTx2tpaJCUloUuXLgDq/o4OHjwIrVZruuBEdEssN0R0g6tXr+KOO+7At99+i99++w3p6en44Ycf8P777xsKwZAhQzB48GCMGzcOu3fvRnp6Onbs2IGdO3cCADp06IDdu3fjyJEjSElJwVNPPYX8/PybvmbHjh3x6KOPYtKkSdi0aRPS09ORmJiIBQsW1DtV2xQKCgpw4sQJ3HvvvTdd5+6770ZhYSHefvvtRp9r6dKl2Lx5M86dO4cZM2agpKTE8DXWs88+C7VajYcffhgnTpzAhQsX8M033yA1NdWk74eI6uPXUkR0A2dnZ0RHR+Pjjz/GpUuXoNVqERgYiGnTpuG1114zrLdx40a8/PLLmDBhAioqKhAWFmY4I+mNN95AWloaRowYAUdHR0yfPh1jx45FaWnpTV931apV+L//+z+89NJLyMnJgaenJ/r169doCWmOn3/+GVFRUfD09LzpOhKJpNHH/7Rw4UIsXLgQp06dQlhYGH766SfDdm3atMGePXswa9YsDBkyBDKZDBERERgwYIDJ3gsR3UgiCIIgdggiouaSSCTYvHkzxo4d2+RtxowZg4EDB+I///lPs183IyMDoaGhOHny5G1P9RASEoIXXnjhtqeiIKI6/FqKiKzehAkTjJrCYODAgZgwYYIZEzXN/Pnz4ezsfNMzyIioeXjkhois2sWLFwEAMpkMoaGhLfa6pjhyU1xcbJjSoW3btlCpVCZMSNR6sdwQERGRTeHXUkRERGRTWG6IiIjIprDcEBERkU1huSEiIiKbwnJDRERENoXlhoiIiGwKyw0RERHZFJYbIiIisiksN0RERGRT/h848ZwBnibFQwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(hm.k_hm, hm.tracer_profile.u(hm.k_hm, m=1e12))\n", "plt.xlabel(\"Scale [h/Mpc]\")\n", "plt.ylabel(\"Normalized Fourier Density\")\n", "plt.xscale(\"log\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can see it's indeed around 0.1 Mpc$^{-1}$h" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally we can see the power spectrum:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAG4CAYAAACpRojiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABqcElEQVR4nO3deVxUVf8H8M/MwLDvq2yCu7iAIeCaG2ZqmpqVLeZSVq4lPyutzNJKy8cyw0cfNVNT09I0y3IjEfcFBPcNUVB2kB1mmJn7+8OaIlwYmOEOM5/368XreebMnTsfvMF8OefccySCIAggIiIiMkNSsQMQERERiYWFEBEREZktFkJERERktlgIERERkdliIURERERmi4UQERERmS0WQkRERGS2WAgRERGR2bIQO4Cx02g0yMjIgIODAyQSidhxiIiIqBYEQUBJSQl8fHwgld6/34eF0ENkZGTA399f7BhERERUB+np6fDz87vv8yyEHsLBwQHA3X9IR0dHkdMQERFRbRQXF8Pf31/7OX4/LIQe4q/hMEdHRxZCREREjczDprVwsjQRERGZLRZCREREZLY4NEZERGRgarUaVVVVYscwKZaWlpDJZPU+DwshIiIiAxEEAVlZWSgsLBQ7iklydnaGt7d3vZa3YSFERERkIH8VQZ6enrC1teV6dHoiCALKy8uRk5MDAGjSpEmdz8VCiIiIyADUarW2CHJzcxM7jsmxsbEBAOTk5MDT07POw2ScLE1ERGQAf80JsrW1FTmJ6frr37Y+869YCBERERkQh8MMRx//tiyEiIiIyGyxECIiIiKzxUKIiIiIdBIYGIjFixfX6xwffvghQkND9ZKnPlgIERERUTXx8fEYMmQIfHx8IJFIsH37drEjGQwLISIiIqqmrKwMISEhWLp0qdhRDI6FEBERUQMrV6ru+1VZpdbrsXUxcOBAfPzxxxg+fPj9v4fycowfPx4ODg4ICAjAihUrqj3/zjvvoFWrVrC1tUWzZs0we/bsB97mrtFoMHfuXPj5+cHKygqhoaHYtWtXnfLrggsqEhERNbDgD3bf97k+rT3w7bgI7eOweftQ8a+C5y+RQa7Y/FpX7eMen+1HQZmy2jE3FgyuZ9p7W7RoEebNm4d3330XW7ZswcSJE9GrVy+0bt0aAODg4IA1a9bAx8cHZ8+exYQJE+Dg4IC33377nuf76quvsGjRIvzvf/9Dp06dsHr1agwdOhTnz59Hy5YtDfI9AOwRIiIiojoYNGgQJk2ahBYtWuCdd96Bu7s79u/fr33+/fffR7du3RAYGIghQ4ZgxowZ+OGHH+57vv/85z945513MGrUKLRu3RqfffYZQkND6z0p+2HYI0RERNTALswdcN/npP9aJDBhdlStjz30Tp/6BdNBx44dtf9fIpHA29tbu/cXAGzevBlLlixBSkoKSktLoVKp4OjoeM9zFRcXIyMjA927d6/W3r17dyQnJxvmG/gTCyEiIqIGZiuv/cevoY6tL0tLy2qPJRIJNBoNAODo0aN44YUX8NFHH2HAgAFwcnLCpk2bsGjRogbLV1sshIiIiEivjhw5gqZNm+K9997Ttt28efO+xzs6OsLHxweHDx9Gr169tO2HDx9GRETEfV+nDyyEiIiIqJrS0lJcu3ZN+zg1NRVJSUlwdXVFQEDAQ1/fsmVLpKWlYdOmTQgPD8fOnTuxbdu2B77mrbfewpw5c9C8eXOEhobi22+/RVJSEjZs2FDv7+dBWAgRERFRNadOnUKfPn/PN4qOjgYAjBkzBmvWrHno64cOHYrp06djypQpUCgUGDx4MGbPno0PP/zwvq+ZNm0aioqK8H//93/IyclBcHAwduzYYdA7xgBAIgiCYNB3aOSKi4vh5OSEoqKi+07yIiIi+rfKykqkpqYiKCgI1tbWYscxSQ/6N67t5zdvnyciIiKzxUKIiIiIzBYLISIiIjJbLISIiIjIbLEQIiIiIrPFQoiIiIjMFgshIiIiMltmsaBiYGAgHB0dIZVK4eLiUm13XCIiIjJfZlEIAXf3PbG3txc7BhERkUmQSCTYtm0bhg0bVudzjB07FoWFhdi+fbvecumKQ2NERERUzfz58xEeHg4HBwd4enpi2LBhuHz5stixDMLoC6H4+HgMGTIEPj4+kEgk96waly5disDAQFhbWyMyMhInTpyo9rxEIkGvXr0QHh5u8M3biIiIGrsDBw5g8uTJOHbsGPbu3Yuqqio89thjKCsrEzua3hl9IVRWVoaQkBAsXbr0ns9v3rwZ0dHRmDNnDhITExESEoIBAwYgJydHe8yhQ4eQkJCAHTt24NNPP8WZM2fu+34KhQLFxcXVvoiIiMzJrl27MHbsWLRr1w4hISFYs2YN0tLSkJCQUO24vLw8DB8+HLa2tmjZsiV27NihfU6tVuPll19GUFAQbGxs0Lp1a3z11VcPfF+FQoFp06bB09MT1tbW6NGjB06ePGmQ7/EvRl8IDRw4EB9//DGGDx9+z+e/+OILTJgwAePGjUNwcDCWL18OW1tbrF69WnuMr68vAKBJkyYYNGgQEhMT7/t+8+fPh5OTk/bL399fv98QERGZvXKlCuVKFf6577lSpUG5UgWFSn3PYzWav4+tUt89trLq4cfqQ1FREQDA1dW1WvtHH32EZ555BmfOnMGgQYPwwgsvoKCgAACg0Wjg5+eHH3/8ERcuXMAHH3yAd999Fz/88MN93+ftt9/G1q1bsXbtWiQmJqJFixYYMGCA9pyGYPSF0IMolUokJCQgKipK2yaVShEVFYWjR48CuNujVFJSAgAoLS3FH3/8gXbt2t33nLNmzUJRUZH2Kz093bDfBBERmZ3gD3Yj+IPdKChTattWxKcg+IPdmPPz+WrHhs3bh+APduN2YYW2bd3Rmwj+YDfe2Vp9hKPHZ/sR/MFuXMst1VtWjUaDN998E927d0f79u2rPTd27Fg899xzaNGiBT799FOUlpZqp6dYWlrio48+QufOnREUFIQXXngB48aNu28hVFZWhmXLlmHhwoUYOHAggoODsXLlStjY2OCbb77R2/fzb436rrG8vDyo1Wp4eXlVa/fy8sKlS5cAANnZ2dreJLVajQkTJiA8PPy+57SysoKVlZXhQhMRETUikydPxrlz53Do0KEaz3Xs2FH7/+3s7ODo6FhtasrSpUuxevVqpKWloaKiAkqlEqGhofd8n5SUFFRVVaF79+7aNktLS0RERODixYv6+4b+pVEXQrXRrFkzJCcnix2DiIhI68LcAQAAG0uZtu3VR5tjfI8gyKSSascmzL476mFt8fexL3Vtiuci/CGVVD/20Dt9ahxbH1OmTMGvv/6K+Ph4+Pn51Xje0tKy2mOJRAKNRgMA2LRpE2bMmIFFixaha9eucHBwwMKFC3H8+HG9ZNOXRl0Iubu7QyaTITs7u1p7dnY2vL29RUpFRET0YLbymh+/cgsp5PeYsXKvYy1lUljKandsXQiCgKlTp2Lbtm2Ii4tDUFCQzuc4fPgwunXrhkmTJmnbUlJS7nt88+bNIZfLcfjwYTRt2hQAUFVVhZMnT+LNN9/U+f1rq1HPEZLL5QgLC0NsbKy2TaPRIDY2Fl27dhUxGRERUeM1efJkrF+/Hhs3boSDgwOysrKQlZWFioqKh7/4Ty1btsSpU6ewe/duXLlyBbNnz37gHWB2dnaYOHEi3nrrLezatQsXLlzAhAkTUF5ejpdfflkf39Y9GX2PUGlpKa5du6Z9nJqaiqSkJLi6uiIgIADR0dEYM2YMOnfujIiICCxevBhlZWUYN26ciKmJiIgar2XLlgEAevfuXa3922+/xdixY2t1jtdeew2nT5/Gs88+C4lEgueeew6TJk3C77//ft/XLFiwABqNBqNHj0ZJSQk6d+6M3bt3w8XFpa7fykNJhH/eu2eE4uLi0KdPnxrtY8aMwZo1awAAMTExWLhwIbKyshAaGoolS5YgMjJSL+9fXFwMJycnFBUVwdHRUS/nJCIi01dZWYnU1FQEBQXB2tpa7Dgm6UH/xrX9/Db6HqHevXvjYbXalClTMGXKlAZKRERERKaiUc8RIiIiIqoPFkJERERktlgIERERkdliIURERERmi4UQERERmS0WQkRERGS2WAgRERGR2WIhRERERGaLhRARERFVIwgCXn31Vbi6ukIikSApKemhr5FIJNi+fTsA4MaNG7V+ndiMfmVpIiIiali7du3CmjVrEBcXh2bNmsHd3V3sSAbDQoiIiIiqSUlJQZMmTdCtWzexoxgch8aIiIhIa+zYsZg6dSrS0tIgkUgQGBiIwMBALF68uNpxoaGh+PDDD0XJqE/sESIiImoggiCgokrd4O9rYymDRCKp1bFfffUVmjdvjhUrVuDkyZOQyWQIDw83cELxsBAiIiJqIBVVagR/sLvB3/fC3AGwldfuI9/JyQkODg6QyWTw9vY2cDLxcWiMiIiIzBZ7hIiIiBqIjaUMF+YOEOV960MqlUIQhGptVVVV9TqnsWAhRERE1EAkEkmth6iMiYeHBzIzM7WPi4uLkZqaKmIi/eHQGBERET1Q37598d133+HgwYM4e/YsxowZA5msfr1MxqLxlaVERETUoGbNmoXU1FQ88cQTcHJywrx580ymR0gi/HvQj6opLi6Gk5MTioqK4OjoaND3UmsEpN8px+XMEpQoqmAnt8Cd8ircKVdi/6Uc3CmvQpC7LQCgVKFCdrECt+6UQy6TwtlWDpVGA5VaQGFFFdQaAfZWFrC2lMFCKoFaI+BOuRLWljK08XaAvbUF7K0skFFYAZVaQOdAF3T0c4anoxVc7eRwsLZAE0drSKXsNCQiqovKykqkpqYiKCgI1tbWYscxSQ/6N67t5zd7hEQS88dV/HImA3klSrjay1FUXoX8MiXUmgfXpSm5pTXaqtRqlCkrarSXKlQoVahqtJ26eafGsWduF9Vok0iAdj6OCHK3R5C7HSQSwN/ZBr1ae8DDgT/URETU+LEQEsmmk+m4dedu8ZJfptS2/7Xcla1chjZNHOFiawlnWzkyCisgCECnAGf4u9rCzsoCEARkl1TC2UaOll4OsJBKYCmTIr9UAaVaAxdbOeQWUqg1AooqqnAjvwwqtQBPByttkXQkJR+ZhRXwcrRGuVKN7JJKZBZWoqJKDUEAzt0uxrnbxdWySwA097RHB18nNPewg4O1JQa294anI4sjIiJqXFgIiWR0l6Y4nX4HNpYW6NbcDW28HeHpaAU3OzksZPUdjnK4Z2v3FjU3zRvXPeiex+aXKnA1uwTFlSqk5pUhNa8Mu89nobCiCoIAXMspxbWcv3un5uw4j3Y+jujSzA2RQa7o6OcEbyeben4fREREhsU5Qg/RkHOEGouc4kqcyyjCudvF2H0+Cxczi3GvET0rCymGhPhgaIgPIoJcYV3PdSyIiBoTzhEyPM4RIlF4Olqjr6M1+rbxwrR+LQEAt++UIyGtEMev5+PQtTzczC+HQqXBloRb2JJwCzaWMrT2dkCAqy3Gdm+KRwJcRf4uiIiIWAiRnvi62MLXxRZDQ3wAABcyi/DbmUzklSqx/3IOsosVSEovRFJ6IXYkZ6CNtwMea+eNAcFeaNvEgXenEZHJ0mg0YkcwWfr4t+XQ2ENwaKz+BEHAxcwSLNxzGSdT81GuVFcbSpNJJejo54RZA9sirKkLZNLa7ZBMRGTMNBoNrl69CplMBg8PD8jl8lrvAE8PJggClEolcnNzoVar0bJlyxp/UNf285uF0EOwENK/wnIlYi/mYPf5LMReyqm2ZICHgxUeb+eNtk0cMKKTL6wb4VL0RER/USqVyMzMRHl5udhRTJKtrS2aNGkCuVxe4zkWQnrCQsiw8korsTL+Oi5llSIx7Q5KKv9e90giAUY+4ocnQ33RpZmrHu6mIyJqeIIgQKVSQa1Wix3FpMhkMlhYWNy3l42FkJ6wEGo4SpUGh1PysPbIDRy4nIt//ofpYmuJVl4O6NbCDS93bwZ7a/YUERHR/bEQ0hMWQuIoV6qw61wWTt4owO7z2Sj4x6KTtnIZBrTzxuPtvfFoSw/YyHlbPhERVcdCSE9YCIlPpdYg7kouPt91CSm5ZdXmFFnKJHCxlePJUB9M6t0CLnY1x4mJiMj8sBDSExZCxkWjEZCYdge/nc3C7vNZuF349x5rMqkE4YEu6NvGE6H+zogIchMxKRERiYmFkJ6wEDJegiDg1zMZWH8sDQVlSlzNqb4hrZ2VDOO7B6F/sBc6+DrxtlUiIjPCQkhPWAg1HukF5dhzIRtrj9xAWkH1W1WbOFnD39UWI8P8MPIRXy7gSERk4lgI6QkLocYpNbcU+y/n4kRqAeKv5qJc+fdtq77ONhjWyQdDQ3zR2vveG9QSEVHjxkJIT1gINX6VVWr8eCodyw9cR2ZRRbVVrd3t5ejbxhMfPNGOt+QTEZkQFkJ6wkLItJQpqhB7KRc7kjJw4EoOqtR3//N3sbXE672a47nIADhaW4qckoiI6ouFkJ6wEDJd2cWVeGfrGRxJyYdSdXfjPnsrC/Rt44lXH22G9r5OIickIqK6YiGkJyyETJ9CpcbPSRlYGX+92p1nL0YG4MOh7bi1BxFRI8RCSE9YCJkPjUbAz0kZeG/7We3k6va+jlgwoiN7h4iIGhkWQnrCQsj8qFQa/O/gdfzvQAqKK1WQSoBerTzw2VMd4eloLXY8IiKqhdp+frPPn+hfLCykmNynBfb9Xy880bEJNAKw/3Iuus7/A7vOZYodj4iI9IiFENF9eDpYI+b5R/D2460hk0igFgRM3JCIL/ZeqbbfGRERNV4shIgeYlLvFjj+Xl8809kPggAsib2K51cew4XMIrGjERFRPbEQIqoFd3trfD4yBIufDYWtXIbjqQUYvOQQVh28LnY0IiKqBxZCRDoY1skX30/oArmFFIIAfPLbRXy17yqHyoiIGikWQkQ6CvF3xrGZ/dCzpTsEAfhy3xWMWnEUKbmlD38xEREZFRZCRHXgai/Hdy9HYtHTIbCTy3Dyxh1ELTqA6ZuSoNFoxI5HRES1xEKIqB6eCvPDrjcfhY+zDQQA25JuY/zaU8gpqRQ7GhER1QILIaJ68ne1xYEZvdGvjScspBLEXc7FgC/j8UtSBnuHiIiMHFeWfgiuLE26uJxVgumbk3AhsxjA3U1cZz8RjGfD/UVORkRkXriy9J8KCwvRuXNnhIaGon379li5cqXYkciEtfZ2wPbJ3TGtbwtIAJQqVHhn6xlM2pCA1LwyseMREdG/mHyPkFqthkKhgK2tLcrKytC+fXucOnUKbm5utXo9e4SorpJvFeLtLWdwNbsEGgGwkEowJMQHk3s3RwsvB7HjERGZNPYI/Ukmk8HW1hYAoFAoIAgCTLz2IyMR4ueM3W8+it/e6Ik+rT2g0gjYdvo2+n8Zjy/2XOZ/h0RERsDoC6H4+HgMGTIEPj4+kEgk2L59e41jli5disDAQFhbWyMyMhInTpyo9nxhYSFCQkLg5+eHt956C+7u7g2Ungho4+2Ib8dFYMXoMFhZSCEAWPLHNUzakIjcEoXY8YiIzJrRF0JlZWUICQnB0qVL7/n85s2bER0djTlz5iAxMREhISEYMGAAcnJytMc4OzsjOTkZqamp2LhxI7KzsxsqPpHWY+28kfB+fzwb7g8LqQS/n8vCY18ewMr467y7jIhIJI1qjpBEIsG2bdswbNgwbVtkZCTCw8MRExMDANBoNPD398fUqVMxc+bMGueYNGkS+vbti5EjR97zPRQKBRSKv/9KLy4uhr+/P+cIkV6dzyjCjB/P4OKfd5d5Olhh82tdEORuL3IyIiLTYBZzhJRKJRISEhAVFaVtk0qliIqKwtGjRwEA2dnZKCkpAQAUFRUhPj4erVu3vu8558+fDycnJ+2Xvz9veyb9a+fjhB1TumNAOy8AQE6JAq+sPYWsIi7ESETUkBp1IZSXlwe1Wg0vL69q7V5eXsjKygIA3Lx5Ez179kRISAh69uyJqVOnokOHDvc956xZs1BUVKT9Sk9PN+j3QObLUibF/0Z3xn9feARejlZIyS3DyOVHcDOft9kTETUUC7EDGFpERASSkpJqfbyVlRWsrKwMF4joXwZ1aIIOvk4Y/c1x3MgvxxNfH8LCkR3xePsmYkcjIjJ5jbpHyN3dHTKZrMbk5+zsbHh7e4uUikh3/q62+OH1rvB1tkFJpQoT1ydi7wVO6iciMrRGXQjJ5XKEhYUhNjZW26bRaBAbG4uuXbuKmIxId54O1vh+QhfYymUQALy56TSOpOSJHYuIyKQZfSFUWlqKpKQk7fBWamoqkpKSkJaWBgCIjo7GypUrsXbtWly8eBETJ05EWVkZxo0bJ2JqoroJcLNF/Ft90K25G8qUaoz99iR+Pn1b7FhERCbL6G+fj4uLQ58+fWq0jxkzBmvWrAEAxMTEYOHChcjKykJoaCiWLFmCyMhIvbw/t9ggMVRWqfHGptPYff7u8FhkkCu+nxAJqdTo/3YhIjIKtf38rlchVFpaWmMhOFMrFlgIkVjUGgFvbDqNX89kAgB6tfLA/0aHwdpSJnIyIiLjZ7B1hFJTUzF48GDY2dnByckJLi4ucHFxgbOzM1xcXOoVmoj+JpNKEPP8I3jqEV9IJcCBK7l4fuUxbstBRKRHOvcIde/eHYIg4I033oCXlxckEkm153v16qXXgGJjjxAZgyPX8vD6+gQUV6rQxMkaHw5thwHteGckEdH9GGxozN7eHgkJCQ9cndmUsBAiY3E9txTj1pzEzfxyAMDLPYIw+4lgkVMRERkngw2NhYeHc7VlIhE087DHuvERcLKxBAB8cygVM7eeQWWVWuRkRESNl849QikpKXj99dfx4osvon379rC0tKz2fMeOHfUaUGzsESJjo1Rp8OGO8/j+ZBoEAWjj7YAvnglBsI+T2NGIiIyGwYbGjh07hueffx43btz4+yQSCQRBgEQigVptWn+dshAiY3Xoah7e3HwaeaVKSACM7toUc59sL3YsIiKjUNvPb533Ghs/fjw6deqE77///p6TpYmoYfRo6Y7fpvXEkK8PIbtEgXVHb8LeygJvDWjNn0siolrSuUfIzs4OycnJaNGihaEyGRX2CJGxq1SqMG7NKRy9ng8AmNa3BaIfM4+bGYiI7sdgk6X79u2L5OTkeoUjIv2xllvg+1e7YM6Qu3eQLfnjGub+cl7kVEREjYPOQ2NDhgzB9OnTcfbsWXTo0KHGZOmhQ4fqLRwR1d647kFQqQV88ttFrD58A9dyS7FuvH62miEiMlU6D409aK8jTpYmEt+Y1Sdw4EouAODdgW3waq/mIiciImp4Bhsa02g09/0ytSKIqDFaOz4CQ0KaAAA+/f0SVh9KFTkREZHx4lbWRCbo6+cewdS+d29omPvrBby9JbnGBslERFSHOUJz58594PMffPBBncMQkf5E928FQQBi9l/DD6duITm9EDun9YSFjH//EBH9Rec5Qp06dar2uKqqCqmpqbCwsEDz5s2RmJio14Bi4xwhauxe/+4Udp3PBgA8FuyFJc91grWlTORURESGZbCVpe/3ZmPHjsXw4cMxevTo+p7OqLAQIlOwIj4F/9lzBUqVBmFNXbDqpc5wsZOLHYuIyGAatBACgLNnz2LIkCHVtt4wBSyEyFScSC3AK2tPorhSBUdrC6wa0xkRQW5ixyIiMgiD3TV2P0VFRSgqKtLX6YhIzyKCXLFlYjfYW8lQXKnCqBXHcDmrWOxYRESi0nmy9JIlS6o9FgQBmZmZ+O677zBw4EC9BSMi/Wvl5YCNE7rg6eVHoVBpMGb1SWx6tQsC3e3EjkZEJAqdh8aCgoKqPZZKpfDw8EDfvn0xa9YsODg46DWg2Dg0RqYop7gSL6w6jqs5pWjiZI3Nr3ZFgJut2LGIiPSmwecImSoWQmSqcksUGLXiKFJyy2BtKcW68RGcM0REJsMgc4SqqqpgYWGBc+fO1TsgEYnLw8EK30/oAju5DJVVGry46gRu3SkXOxYRUYPSqRCytLREQEAAt9IgMhGejtb4/tUusLKQQqnW4PmVx5FZVCF2LCKiBqPzXWPvvfce3n33XRQUFBgiDxE1sI5+ztg/oxcCXG2RVlCOUSuO4TZ7hojITNRpZelr166hqqoKTZs2hZ1d9btNuLI0UeN0u7ACo1YcRXpBBeQyKTZOiETnQFexYxER1UltP791vn3+ySefhEQiqVc4IjI+vs422DShC/ouOgCFSoMXvzmOvdN7wd+Vd5MRkeniXWMPwR4hMjdn0gsxauUxlCvV8HW2wcYJkWjqxnWGiKhxMdjK0s2aNUN+fn6N9sLCQjRr1kzX0xGRkeno74z9M3qjmYfdn8Nlx3Alq0TsWEREBqFzIXTjxo173jWmUChw69YtvYQiInF5OVpj06td0MLTHplFlXj8q3jEX8kVOxYRkd7Veo7Qjh07tP9/9+7dcHJy0j5Wq9WIjY2tseo0ETVeng7W2PhKJHp+vh8KlQavr0/A3uhe8HW2ETsaEZHe1HqOkFR6t/NIIpHg3y+xtLREYGAgFi1ahCeeeEL/KUXEOUJk7q5kl+Dp5UdRVFGFZu522PxaV3g4WIkdi4jogfQ+R0ij0UCj0SAgIAA5OTnaxxqNBgqFApcvXza5IoiI7m7U+vsbPeHrbIPreWV4afUJ5JUqxI5FRKQXOs8RSk1Nhbu7uyGyEJGR8nG2wfpXIuFub4WLmcXotXA/cksqxY5FRFRvOhdC06ZNw5IlS2q0x8TE4M0339RHJiIyQkHudlgyKhQAUKZQ4/mVx1FZxe12iKhx07kQ2rp1K7p3716jvVu3btiyZYteQhGRcerWwh0LRnSApUyCqzmlmPb9aajUGrFjERHVmc6FUH5+frU7xv7i6OiIvLw8vYQiIuM1KiIAa8dFQG4hxZ4L2ZjxYzJUKhZDRNQ46VwItWjRArt27arR/vvvv3NBRSIz0a2FO/77/COwkEqwPSkDA5cchEbDYoiIGh+d9xqLjo7GlClTkJubi759+wIAYmNjsWjRIixevFjf+YjISEUFe2FS7+ZY8sc1XM0pxYwfz2DRMyHci5CIGhWdC6Hx48dDoVDgk08+wbx58wAAgYGBWLZsGV566SW9ByQi4xX9WGtkFFZiS+It/HT6Npq62eGNqJZixyIiqrV6bbqam5sLGxsb2Nvb6zOTUeGCikQP982hVMz79QIA4N1BbfDqo81FTkRE5s5gm64CgEqlwr59+/DTTz9pV5nOyMhAaWlp3dISUaP2co8gRPdvBQD49LdLiN6cJG4gIqJa0rkQunnzJjp06IAnn3wSkydPRm7u3Y0YP/vsM8yYMUPvAYmocZjatwW6NnMDAPx0+ja2JKSLnIiI6OF0LoTeeOMNdO7cGXfu3IGNzd+bLw4fPhyxsbF6DUdEjYdEIsGGVyIQ3ORuF/Ssn87i4FXuWE9Exk3nQujgwYN4//33IZfLq7UHBgbi9u3begtGRI2PVCrFL1N7YFAHb1SpBbz2XQJOp90ROxYR0X3pXAhpNBqo1TWX1b916xYcHBz0EoqIGi+ZVIIvnw1FjxbuKFeqMXLZUcRezBY7FhHRPelcCD322GPV1guSSCQoLS3FnDlzMGjQIH1mI6JGyspChv+NDoODtQXUgoCJGxJx60652LGIiGrQ+fb5W7duYcCAARAEAVevXkXnzp1x9epVuLu7Iz4+Hp6enobKKgrePk9Udym5pXgy5jBKFSo0c7fDj693hZu9ldixiMgM1Pbzu07rCKlUKmzatAlnzpxBaWkpHnnkEbzwwgvVJk+bChZCRPWTWVSBkcuO4nZhBUL8nLBxQhfYWem8lisRkU4MWgiZExZCRPWXkluKkcuO4E55FTzsrXDg7d6wlbMYIiLDMeiCipcvX8aUKVPQr18/9OvXD1OmTMGlS5fqHJaITFtzD3ssea4TACC3VIFn/3cMGg3/BiMi8elcCG3duhXt27dHQkICQkJCEBISgsTERHTo0AFbt241REYiMgE9W3rgzX4tIQFw9nYRPtvFP56ISHw6D401b94cL7zwAubOnVutfc6cOVi/fj1SUlL0GlBsHBoj0q8tCbcw48dkAMD7g9vilZ7NRE5ERKbIYENjmZmZ99xl/sUXX0RmZqaup2sQw4cPh4uLC0aOHCl2FCKzNzLMD28/3hoA8PHOi/hwxzmRExGROdO5EOrduzcOHjxYo/3QoUPo2bOnXkLp2xtvvIF169aJHYOI/jSxV3P0a3N3qY01R27ix1Pcl4yIxKHzbRtDhw7FO++8g4SEBHTp0gUAcOzYMfz444/46KOPsGPHjmrHGoPevXsjLi5O7BhE9CeJRILlL4bh0YX7kVlUiY9+uYB2Pk4I9uHwMxE1LJ3nCEmltetEkkgk99yKQ1fx8fFYuHAhEhISkJmZiW3btmHYsGHVjlm6dCkWLlyIrKwshISE4Ouvv0ZERES1Y+Li4hATE4MtW7bo9P6cI0RkOBVKFcatOYlj1wvg6WCFnyZ1g5+LrdixiMgEGGyOkEajqdWXPoogACgrK0NISAiWLl16z+c3b96M6OhozJkzB4mJiQgJCcGAAQOQk5Ojl/cnIsOxkVvgf6M7o7WXA3JKFBgacxhpBdyKg4gaTp3WEWpIAwcOxMcff4zhw4ff8/kvvvgCEyZMwLhx4xAcHIzly5fD1tYWq1evrtP7KRQKFBcXV/siIsNxsrHEmvHhsJXLUFCmxOAlB1GhVIkdi4jMRK0LoaNHj+LXX3+t1rZu3ToEBQXB09MTr776KhQKhd4DPohSqURCQgKioqK0bVKpFFFRUTh69Gidzjl//nw4OTlpv/z9/fUVl4juo4mTDRY+3RESACWVKkzfnAw1F1wkogZQ60Jo7ty5OH/+vPbx2bNn8fLLLyMqKgozZ87EL7/8gvnz5xsk5P3k5eVBrVbDy8urWruXlxeysrK0j6OiovD000/jt99+g5+f3wOLpFmzZqGoqEj7lZ7Ou1mIGsLgDj5YMy4ccpkUu85nYe4v58EdgIjI0Gp911hSUhLmzZunfbxp0yZERkZi5cqVAAB/f3/MmTMHH374od5D1te+fftqfayVlRWsrLg7NpEYerX2xBfPhmDKxtNYe/QmyhQq/OeZULFjEZEJq3WP0J07d6r1vBw4cAADBw7UPg4PD2/w3hN3d3fIZDJkZ2dXa8/Ozoa3t3eDZiEi/Xiiow9mPNYKALAl8TY+3XlB5EREZMpqXQh5eXkhNTUVwN25OYmJidp1hACgpKQElpaW+k/4AHK5HGFhYYiNjdW2aTQaxMbGomvXrg2ahYj0Z1Lv5mjlZQ8AWH34Bk6kFoiciIhMVa0LoUGDBmHmzJk4ePAgZs2aBVtb22orSZ85cwbNmzfXe8DS0lIkJSUhKSkJAJCamoqkpCSkpaUBAKKjo7Fy5UqsXbsWFy9exMSJE1FWVoZx48bpPQsRNQypVIpfpvRAr1YeUGkETFh3CtdySsWORUQmqNYLKubl5WHEiBE4dOgQ7O3tsXbt2mq3tPfr1w9dunTBJ598oteAcXFx6NOnT432MWPGYM2aNQCAmJgY7YKKoaGhWLJkCSIjI/Xy/lxQkUg8FUo1nl91DKfTCtHEyRrfvRyBFp4OYsciokagtp/fOq8sXVRUBHt7e8hksmrtBQUFsLe3h1wur1tiI8VCiEhc+aUKPPH1IWQWVcJWLsPhmX3gYssbGojowQy2srSTk1ONIggAXF1dTa4IIiLxudlb4ZPhHSABUK5UY+rGJKjUGrFjEZGJMPqVpYmI+rbxxIKnOsDKQopD1/Lw0S8XuMYQEekFCyEiahSeDQ/A4mdDIZEA3x27iZg/rokdiYhMAAshImo0BnZogvcGtQUALNp7BV/uvSJyIiJq7HQqhKqqqjB+/HjtekJERA3t5R5BCHK3BQB8/cdVJKbdETkRETVmOhVClpaW2Lp1q6GyEBE9lEQiwa9Te6CVlz00AvDK2lO4mV8mdiwiaqR0HhobNmwYtm/fboAoRES1Y2dliW2TuqO9ryMKypQY9+1J5JcqxI5FRI1QrTdd/UvLli0xd+5cHD58GGFhYbCzs6v2/LRp0/QWjojofuysLLB6TDiGLT2M63lliPriAI7M7Asbuc6/1ojIjOm8oGJQUND9TyaR4Pr16/UOZUy4oCKRcdt7IQsT1iUAALo0c8X3E7pAIpGInIqIxFbbz2+d/3TiRGkiMib9g73xeq9mWBF/HceuF+DLvVcQ/VhrsWMRUSNR59vnlUolLl++DJVKpc88REQ6mzmwLeaP6AAAWPLHNfx4Kl3kRETUWOhcCJWXl+Pll1+Gra0t2rVrp90FfurUqViwYIHeAxIR1caz4QGY1Ls5AODtrWewJYHFEBE9nM6F0KxZs5CcnIy4uDhYW1tr26OiorB582a9hiMi0sWMx1oj0M0WggC8veUMLmUWix2JiIyczoXQ9u3bERMTgx49elSbkNiuXTukpKToNRwRkS6kUgk2TIiEjaXs7hpD604ht4S31RPR/elcCOXm5sLT07NGe1lZGe/UICLR+Trb4uDbfRDoZotbdyrwyrpTqFCqxY5FREZK50Koc+fO2Llzp/bxX8XPqlWr0LVrV/0lIyKqI3cHK3w7LgLOtpZITi/EqBVHodFwt3oiqknn2+c//fRTDBw4EBcuXIBKpcJXX32FCxcu4MiRIzhw4IAhMhIR6SzI3Q5fPh2KcWtPIvlWEV5bn4CVL3UWOxYRGRmde4R69OiBpKQkqFQqdOjQAXv27IGnpyeOHj2KsLAwQ2QkIqqTPm090aOFOwBg74VsrD92U+RERGRsdF5Z2txwZWmixk2l1uDLvVewNC4FMqkE34zpjN6ta85zJCLTUtvPb517hF566SV8++23JreVBhGZJguZFDMGtMbIMD+oNQImrk/E6bQ7YsciIiOhcyEkl8sxf/58tGjRAv7+/njxxRexatUqXL161RD5iIjqTSKR4NPhHdDOxxEVVWo8t+IYMu5UiB2LiIyAzoXQqlWrcOXKFaSnp+Pzzz+Hvb09Fi1ahDZt2sDPz88QGYmI6k1uIcX/9W8FAKhUaTBuzUmUKrhFEJG5q/NeYy4uLnBzc4OLiwucnZ1hYWEBDw8PfWYjItKrvm298PVzoXC3l+Nydgkmb0hElVojdiwiEpHOhdC7776Lbt26wc3NDTNnzkRlZSVmzpyJrKwsnD592hAZiYj0ZkiIL1aPDYeNpQwHruTi3W1nwXtGiMyXzneNSaVSeHh4YPr06RgxYgRatWplqGxGgXeNEZmmvRey8eq6UxAATOrdHG8/3kbsSESkRwa7a+z06dN47733cOLECXTv3h2+vr54/vnnsWLFCly5cqVeoYmIGkrPlu5wt5cDAP4bl4IdyRkiJyIiMdR7HaHk5GR8+eWX2LBhAzQaDdRq09rThz1CRKYrLb8cM35MxokbBZDLpFj/SiQiglzFjkVEelDbz2+dt9gQBAGnT59GXFwc4uLicOjQIRQXF6Njx47o1atXvUITETWkADdbfP9qF0zakIDd57MxYd0pbHm9K1p6OYgdjYgaiM49Qi4uLigtLUVISAh69eqF3r17o2fPnnB2djZQRHGxR4jI9FUo1Xh+1TGcTiuEi60lfn/jUXg7WYsdi4jqwWA9QuvXr0fPnj1ZFBCRybCRyzCovTdOpxXiTnkVRn9zHFsndYOjtaXY0YjIwHSeLD148GBtEXTr1i3cunVL76GIiBrauO5B6NnSHQ7WFriaU4rXv0uAUsU1hohMnc6FkEajwdy5c+Hk5ISmTZuiadOmcHZ2xrx586DR8JcGETVOFjIp1o2PwPcTusBOLsORlHy8tSUZGg3XGCIyZToXQu+99x5iYmKwYMECnD59GqdPn8ann36Kr7/+GrNnzzZERiKiBiGRSNDe1wnLXgyDhVSCn5MysGDXJbFjEZEB6TxZ2sfHB8uXL8fQoUOrtf/888+YNGkSbt++rdeAYuNkaSLzU6XWoNfn+5FRVAkAeG9QW0x4tJnIqYhIFwZbULGgoABt2tRcgbVNmzYoKCjQ9XREREbHUibF+B5BsJXLAACf/HYRP5xMFzkVERmCzoVQSEgIYmJiarTHxMQgJCREL6GIiMT2co8gHH6nL17rdbcnaOZPZ/D72UyRUxGRvul8+/znn3+OwYMHY9++fejatSsA4OjRo0hPT8dvv/2m94BERGKQSCRwsZNj5uNtUFRehU0n0zFt02mstrZAz5YeYscjIj3RuUeoV69euHLlCkaMGIHCwkIUFhZixIgRuHz5Mnr27GmIjEREopFIJBgZ5gcrCymq1AJe+y4BiWl3xI5FRHqiU4/QjRs3sHfvXiiVSowaNQrt27c3VC4iIqNhZSGDIAiwlctQrlRj3Lcnsfm1LmjjzRsoiBq7Wt81tn//fjzxxBOoqKgAAFhYWGD16tV48cUXDRpQbLxrjIgA4Nj1fDTzsMPr3yUgMa0Q7vZW+H5CJPclIzJSer9rbPbs2ejfvz9u376N/Px8TJgwAW+//bZewhIRGbsuzdzg6WCNb8dGILiJI/JKFRi14hguZ5WIHY2I6qHWPULOzs44cuQIgoODAQDl5eVwdHREdnY23NzcDBpSTOwRIqJ/++FkOt7ffg5KtQaudnKsfzkSwT78/UBkTPTeI1RcXAx3d3ftY1tbW9jY2KCoqKh+SYmIGhGlSoMVB69DqdbAw8EKBWVKPL/qGM7d5u9CosZIp8nSu3fvhpOTk/axRqNBbGwszp07p23794rTRESmRG4hxYrRYfgp8TbGdw/E+LWnkJReiBdWHcf6lyPRwc/p4SchIqNR66ExqfThnUcSiQRqtbreoYwJh8aI6EGKK6swdvUJJKYVwsHaAmvHR+CRABexYxGZPb0PjWk0mod+mVoRRET0MA5WFujV2gNtvB1QUqnCcyuO4TeuQE3UaOi8oCIREf1t/bGb+HLvVeQUV6JnS3coVBpM2pCIZXEp0HFPayISAQshIqJ6GBnmjxB/Z0Q/1hrfjg3H2G6BAIDPdl3CzK1noVRpxA1IRA9U6zlC5opzhIjoYdQaATKpRPt4zeFUzP31AjQC0K25G5a9EAYnW0sRExKZH73PESIionv7ZxGkVGnQzMMeq8Z0hp1chiMp+Rix7DAuZRWLmJCI7oeFEBGRnpQrVXhu5TGM/fYELKRS/Ph6NzRxskZKbhmGfH0IS/dfg0rNoTIiY1LnQighIQHr16/H+vXrkZiYqM9MRESNko2lDM097GBnZQGNICDYxxE7pvRAVFsvVKkFLNx9GU8tP4prOaViRyWiP+k8RygnJwejRo1CXFwcnJ2dAQCFhYXo06cPNm3aBA8PD0PkFA3nCBGRLhQqNbKKKtHUzU7bJggCfkq8jQ9/OY+SShWsLKR4a0BrjOseVG1YjYj0x2BzhKZOnYqSkhKcP38eBQUFKCgowLlz51BcXIxp06bVK7ShDB8+HC4uLhg5cqTYUYjIxFlZyKoVQSWVVQCAp8L8sGf6o3i0lQcUKg0+3nkRTy07gv2Xc3ibPZGIdO4RcnJywr59+xAeHl6t/cSJE3jsscdQWFioz3x6ERcXh5KSEqxduxZbtmzR6bXsESKiurqSXYIJ607h6TA/TOnbEsDd3qFNJ9Px8a8XUKa8uwhtcBNHTOzdHIM6NGEPEZGeGKxHSKPRwNKy5m2glpaW0GiMcxJg79694eDgIHYMIjIziTfv4GZ+OTafSke5UgXg7lZEz0UE4I8ZvTGhZxBs5TJcyCzG1O9PI+qLA9h0Ig2VVVyln6ih6FwI9e3bF2+88QYyMjK0bbdv38b06dPRr18/nQPEx8djyJAh8PHxgUQiwfbt22scs3TpUgQGBsLa2hqRkZE4ceKEzu9DRNTQRkUE4KOh7fDz5B6wlVff49rL0RrvDQ7GkZl9MT2qFZxtLZGaV4aZP51FyEd78NLqE1h18DquZpdw6IzIgHTafR4AYmJiMHToUAQGBsLf3x8AkJ6ejvbt22P9+vU6BygrK0NISAjGjx+PESNG1Hh+8+bNiI6OxvLlyxEZGYnFixdjwIABuHz5Mjw9PQEAoaGhUKlUNV67Z88e+Pj46JRHoVBAoVBoHxcXc+0PIqq7MX+uNH0/zrZyvBHVEq/0DML3J9Lw7eEbuF1YgfgruYi/kouPd15EEydr9GjhjkB3O/g626CJkzV8nG3g7WQNSxlXQSGqjzqtLC0IAvbt24dLly4BANq2bYuoqKj6h5FIsG3bNgwbNkzbFhkZifDwcMTExAC4OzTn7++PqVOnYubMmbU+d1xcHGJiYh46R+jDDz/ERx99VKOdc4SIqL4OXc3DL8kZmD+iA6T3mQskCAKu5ZTiwJVcHLiSixOpBVDcZ5sOiQRwsrGEpUwKS6kElhZSWEglsJRJ0aWZG6b3bwUnG65oTeaptnOEdO4RWrduHZ599ln0798f/fv317YrlUps2rQJL730Ut0S34NSqURCQgJmzZqlbZNKpYiKisLRo0f19j7/NGvWLERHR2sfFxcXa3u+iIjqKq9UgVfWnURllQadApwxKiLgnsdJJBK09HJASy8HvNKzGSqr1DieWoCEGwW4VViBzMJKZBZVIKOoEkqVBoXlVfc8z6WsEvx2NhNzn2yHx9s3MeS3RtSo6VwIjRs3Do8//rh2WOovJSUlGDdunF4Loby8PKjVanh5eVVr9/Ly0vZG1UZUVBSSk5NRVlYGPz8//Pjjj+jates9j7WysoKVlVW9chMR/Zu7vRXmDGmHhJt3MKyTb61fZ20pQ69WHujVqvoabSq1Bolpd6ARBDhay6HSaHDoai4+330FACABkFOiwOvrEzGgnRfmPtkeXo7W+vyWiEyCzoWQIAiQSGp26d66dQtOTk56CaVv+/btEzsCERGeiwjAqHD/e/4OfRCNRsD1vDKcvV2IM7eKcO52Ec5nFKNcqcbsJ4Lxco8gAIClTIrPd1+Bg5UFShR/z5vcfT4bR67lY+agNnguPOC+w3JE5qjWhVCnTp0gkUggkUjQr18/WFj8/VK1Wo3U1FQ8/vjjeg3n7u4OmUyG7Ozsau3Z2dnw9vbW63sRETWEfxZBv57JQJ/WnrCz+vv3qSAIUKg0sLaUAQAuZRVj5LKjKFXUvCHEVi6rdqt9S097JM7uD2cbS+xIzsCC3y8hq7gSAFCiUOG9beewLfE2Ph3RAa28uKQIEaBDIfTXBOakpCQMGDAA9vb22ufkcjkCAwPx1FNP6TWcXC5HWFgYYmNjte+v0WgQGxuLKVOm6PW9iIga0hd7LmPJH9fQr60nngv3x5lbRUi6VYQztwoxLNQXHw5tBwAIcLVFuVIFa0sp2vk4oYOvEzr63f3fZh721RZgtJBJ4WonBwAM6+SLx9p5YfmB61gedw2+LrbILq7EqZt3MOirg3j10WaY2rclbOQyUb5/ImNR60Jozpw5AIDAwEA8++yzsLbWz1hzaWkprl27pn2cmpqKpKQkuLq6IiAgANHR0RgzZgw6d+6MiIgILF68GGVlZRg3bpxe3p+IqKFVVqlx6sYdAEDsxRzEXsyp9vz5jCLt/7eVW2BvdC80dbWFhY63ytvKLRDdvxWeDfeHSq2BhUyKD3ecx94L2fhvXAp2JGfg42Ht0bu158NPRmSi6nT7vD7FxcWhT58+NdrHjBmDNWvWALi7dtHChQuRlZWF0NBQLFmyBJGRkQ2Sj1tsEFFdCIKAW3cqkJh2B4k378BSJsX7TwRrn+/yaSyyiishkQDNPewR4ueMUH8ndPRzRpsmDrCyMFxPzbhvT2D/5Vzt4yc6NsHsJ4I5mZpMSm0/v3UuhKRS6QMn+qnVprU0PAshIqqt02l3cPJGARJvFiIh7Q5yS/5enNXVTo6E96O0vz9/P5sJRxtLdPBzgqO1JTQaAYUVVdqhLUPadvoWPtl5EXmlSm2bjaUMbz/eGqO7NNW554nIGBlsHaGffvqpWiFUVVWF06dPY+3atfdciJCIyBQVVVThYmYxujRz07Z9+ttFnPxzyAsALKQStPN1wiMBzugU4AKNAMj+/PU5sEOTaueavjkJWUWV2Dqxm8Hn7Qzv5IfHgr2xLC4Fy+NToFILqKhS46NfLuD7E2lY8FRHPBLgYtAMRMZCb0NjGzduxObNm/Hzzz/r43RGgz1CRAQABWVKnEjNx7HrBTh2PR+Xs0sAAEkfPKZdvXlZXAoS0+4grKkLwpq6oIOvk/burwfJKKzAkK8PoVShwrrxEYj8R3FlaOkF5Zj/+0X8djarWvtzEf54e0AbuDRADxWRIRhsaOx+rl+/jo4dO6K0tFQfpzMaLISIzNuWhFtYEZ+CK9k1f7cFutli+egwtPGu/++GkzcKYGMpQ3tfcdZjO3mjAJ/vugRvR2v8ciYTAOBsa4l3Hm+DZzv7c+0hanQMNjR2LxUVFViyZAl8fWu/WioRkTEpV6pw8sYdHLmWh+ciAhDobgcAUKo02iKolZc9IoPcENnMFRFBrvB00N/k4vBAV72dq67v/+Pr3QAAo7sWYPb2c7icXYJZP53FmiM38J+RIejgZ5yL5hLVh86FkIuLS7U5QoIgoKSkBLa2tnXafZ6ISAyCIOBCZjHir+Qh/kouTt0sQJX6bge5j7ONthDq19YTy2wfQUSQK9zsG2b7nVt3yvH+9nOYP6IDmjjZNMh7/lNEkCuWPt8JA746CLVGwOWsEgyJOYSRYX54f3BbONtyuIxMh85DY2vXrq32WCqVwsPDA5GRkXBxMb3JdRwaIzI9524XYdyak9Xu6gIAX2cbdGvuhqfC/KpNgm5oz688hiMp+ejXxhPfjA0XLUd6QTk+3HEesZf+XufIVi7D+4PbYhS36iAj1+BzhADg3LlzaN++vb5OZxRYCBE1btdzSxF7MQeudnI8FeYHAChVqNBp7h5YSKXo2twNj7Z0x6OtPBDkbqfzPmCGcDO/DO9vP4cFT3WEr3PD9wj92/Hr+Xhn6xncyC/XtrVt4oDPnuqIjn7O4gUjeoAGK4RKSkrw/fffY9WqVUhISOA6QkQkKo1GQNKtQuw6l4W9F7KRmlcGAGjn44id03pqjzt7qwitvO0NunChKVFrBPxwMh3zdl5AhVINAYBEAowKD8DbA1rz7jIyOgYvhOLj4/HNN99g69at8PHxwYgRI/DUU08hPFy8blxDYCFE1Hgs3H0JWxNuazcaBQBLmQRdmrmhXxtPjOkWaBQ9Pro6c6sQUolEtDvK/qlUocKZ9EL8cCod25MyAADWllJE92+Fl3s0q7b3GZGYDHLXWFZWFtasWYNvvvkGxcXFeOaZZ6BQKLB9+3YEBwc//ARERHoiCALO3S6udidTekEFsoorYW9lgb5tPDGgnTcebeUOB2tLEZPWz8GruXh57Sl42Fvhl6k9GmTl6Qext7JAtxbu6NbCHc9FBGDq96eRU6LAp79dwurDN/DVs6ENug4SUX3VukdoyJAhiI+Px+DBg/HCCy/g8ccfh0wmg6WlJZKTk022EGKPEJFxuZJdgu2nb2NHcgZu3anA7jcfRWtvBwBAcnoh8koV6N7CvVYLGTYGxZVVeDLmMJp72GHRM6HaxRuNxdGUfLy5+TSyi/+eeN6tuRsWjwrV6/ICRLrS+9CYhYUFpk2bhokTJ6Jly5badhZCRGRod8qU2Hb6Nn5MuIWLmcXadnsrC3z2VEcM7tjkAa9u/HJKKuFuZ2W0d2mpNQJWH0rFf/ZchkKlAQDIpBJM6BmE/3usNSy5dxmJoLaf37X+r/PQoUMoKSlBWFgYIiMjERMTg7y8PL2EJSK6n8S0O4j8NBZzf72Ai5nFsJRJ0D/YCzHPd8Kp96NMvggCAE8H62pFUH6p4gFHNzyZVIIJjzbDyfejMCzUBxLcLY6WH7iOgV8dxKGr/Kwg46XzZOmysjJs3rwZq1evxokTJ6BWq/HFF19g/PjxcHBwMFRO0bBHiKhh5ZYokH6nXLvpp1KlQZf5sWjiZI1R4f4YEuJjtgv6aTQCFsdexbeHUrFtcne08LQXO9I9peSWYurG07hdWIGiiioAQP9gL7w3qA0C3Y0zM5meBrl9/vLly/jmm2/w3XffobCwEP3798eOHTvqejqjxEKIqGGcvVWEb4+k4tfkTHg5WSFuRh/tHUg5JZWcbwJApdbghVXHcTy1AO883gYTezcXO9IDFVVU4cu9V/DdsZtQawRIAAzv5ItPR3QwmTlcZLwadEFFtVqNX375BatXr2YhRES1VqXWYNe5LKw5cgMJN+9o20P9nfG/0WHwcmTx82+5JQocScnDk6GNZ2/H5PQ7GLn8qHYLE2tLKd4d2BajuzZtlMsZUOMgysrSpoiFEJFh7L+Ug/e3n8PtwgoAd9f7GdyhCcZ0C0SnANPbrsfc3SlTYPrmZMRdydW2+ThbI+a5TnikqbgbzpJpatDd54mIdOVgbYHbhRVwt5fjhcimeCEyAJ7sAdJJZZUaH/x8DsM7+aFrc+Neu8fFzgprxkfg7K1CTNqYiPSCCmQUVmLEsqMY2N4bC0Z0hJOtcS0NQOaBPUIPwR4hovorqqjCt4dTAQBvRrXStv92NhN923hyvkgdfbn3Cr6KvQovRysceKtPo/l3FAQBP5xKx0e/XEC58u62TK52csx4rDWeDffn6tSkFxwa0xMWQkR1V1mlxvpjNxGz/xoKy6tgYynD4Zl9RV8d2VRUKNV49btTmNirObq1cBc7js6UKg02nUzDd0dv4mpOKQDAy9EKHwwJxuAOPiKno8aOhZCesBAi0p1aI2Db6dv4cu8V7Ryglp72eCOqJQa2b8K/+KmaKrUG3x29iYV7LqPizx6iVl72WP5iGJp58HZ7qhsWQnrCQohIN+duF+H/fkjG5ewSAEATJ2tMj2qFEY/4woIrDBtUQZkS5zOK0LOlh9hR6uT49Xzt3mUAIAEwqIM3Fj0dAms5p7SSbvS+sjQRUW242slxI78MjtYWmDWwDfbP6I1nwv1ZBBnY7cIKPLHkICasO1VtG5LGJLKZG46/2w+zBraBXCaFAGDn2Sx0/Ggvlh9IAf9uJ0PgbyYiqheFSo29F7K1j32cbbD8xTAcfLsvXuvVvNFM4G3svB2t0cLLAT5ONo166FEikeC1Xs2RPKc/BndoAgkApVqDBb9fwgurjuPKnz2NRPrCobGH4NAY0f0dupqHD34+h+t5Zfjhta6ICOJ6MGIqKq+CRAo4WpvObeipeaV4fX0iUnJKodIIkEkleDEyAJP6tOCCm/RAXEeIiAymqLwKc3acw/akDACAu70VShVVIqeif6/DU1JZBYdGXhQFudtj95uPIi2/HB/vvIA9F7Kx9uhNfHfsJkaFB2Duk+047Er1wv96iEgnB6/mYsDieGxPyoBUAoztFog/ZvRC3zZeYkejf9h7IRs9P99vMju/B7jZYsVLnbFufATkFlJoBGDjiTSEfLQHW06lix2PGjEWQkRUa4v2XMbob04gq7gSzdzt8NOk7vhwaDuTGooxFfsuZKOwvAprj94QO4pePdrKAyff7YeeLe+um1SmVGPGljPotXA/LmQUiZyOGiMOjRFRrTX/c02Xl7o2xayBbWEj50RoY/XRk+3QzMMO43sEiR1F75xs5fju5chq23XczC/HoCWH8ETHJvjP0yGcpE+1xsnSD8HJ0mTOBEFAdrEC3k7W2sfnM4rR3tdJ5GREf9t44ibm/XIRFVV3F2P0dbbBu4PaYlAHb+5ub8a4oKKesBAic1WmUOHtLWeQcPMOfpnaAx4OVmJHojoSBAHfHbsJR2tLDOvkK3Ycg1CpNfj6j2v48VQ6MooqAQDNPOzw7sC2iArm/DVzxAUViajOUvPKMPy/h7HzbCbyyxQ4nXZH7EhUD7vOZeGDn89j5k9nkF5QLnYcg7CQSTG9fyvE/l9vvNGvJeQWUlzPLcMr605h0FfxSMsvEzsiGSn2CD0Ee4TI3PxxKRtvbEpCSaUKng5WWPZiGMKauogdi+pBrRHw6rpT6NrcDS/3CDKL4aLj1/MxZeNp5Jb+vV3HkBAfLBzZEVacP2QWODSmJyyEyFxoNAK+/uMaFsdegSAAnZu64L8vPAJPLlpnEgRBMIsC6J8EQcDyAyn4cu9VKNUaAICVhRTvPN4a43s0EzkdGRqHxohIJysOXseX++4WQS91bYqNE7qwCDIh/yyCVGoNjl3PFzFNw5BIJJjYuwWS5zyGge29IQGgUGkw99eLGPftCdzkcBmBPUIPxR4hMhdFFVV4cdVxvNglAM+GB4gdhwykQqnGuDUncCK1ABsndEGXZm5iR2ow13NL8fr6BFzNKYUgAHKZFK/0DMKEns3gYicXOx7pGXuEiOihKv+83RgAnGwssX1ydxZBJs7aUgofZxvYWMpQVGFe26I087DHnum9sOfNR9GzpTuUag3+G5eCsI/34oOfz0H95/AZmRf2CD0Ee4TIVGUVVWLM6hN4urMfXunJ+RLmpEKpRkZRhXaBTHMkCAL2XczB1O8TUVl1twCys5Jh3pPtMeIRP5HTkT6wR4iI7ut6bimeWnYEl7NLsOpgKkoqzatnwNzZyGXViqAqM+wJkUgk6B/sheOzotC9xd3hwTKFGtE/JKPPf+JwKbNY5ITUUFgIEZmZlNxSPPO/o7hdWIEgdzv8+HrXRr9DOdVdSm4pBi85iN/OZoodRRROtpbY8EoX/DylO/xcbADcXUfr8a8OYtKGRChU6oecgRo7FkJEZuTWnXK8uOo48kqVCG7iiB9f7wp/V1uxY5GItiXexpXsUvxn92WozLBn6C8hfs449E5ffDKsPWws7340/nY2E499GY+9F7LBWSSmi3OEHoJzhMhU5JRU4pnlR3EjvxzNPezww2td4WbPbTPMnUqtwee7L+OVnkHwdOByCQBQpdLgk98uYufZTOSW3F2Qsb2PI6Ifa4W+bbhdR2PBBRX1hIUQmYoNx2/ivW3n4Otsgy0Tu6KJk43YkYiMWqlChaX7r2HVweuoUt/9qGzv44jlo8Pg58KeVGPHQkhPWAiRKdlw/Ca6N3dHoLud2FHISJ26UYC0gnLeOfUPJ2/kY+L6ROSVKgHc3a5jWCcffP5UCCwtOMPEWLEQ0hMWQtSYVVapoREE2MotxI5CjcCZW4UY/t8jkEkk+GlSN7T3dRI7klGJ2X8VX+27qu0dsraUYtbANhjTLUjkZHQvvH2eyMxpNAKmb07C6G9OmN3CeVQ3HXyd0K+NJwa090YQew1rmNKnJZI/eAwDgu/OE6qs0mDOjgsY9+0J3LpTLnI6qiv2CD0Ee4Sosfo69ioW7b0CS5kEm17tgrCmrmJHokagskoNKwup2W3QqqtrOSV4fX0CUnLKIODuZq6v9WqO13o2g501e2CNAYfG9ISFEDVG+y5kY8J3pyAIwIIRHTAqgttmUN1czy1FMzNegfphzmcUYe4vF3A8tQAAIJNK8FKXppj9RFtIpRx0EROHxojM1LWcUry5OQmCAIzu0pRFENWJRiNg3q8XEPXFARy+lid2HKPVzscJm17tgv++8Ahs5TKoNQK+PXIDIR/txY6kDLHjUS2wECIyIcWVVXh13SmUKlSICHTF7CeCxY5EjZRUKkFJZRU0ApCUXih2HKMmkUgwqEMTHH6nL7o0uzsEXaJQYdqm0+i3KA7XckpFTkgPwqGxh+DQGDUm074/jR3JGWjiZI1fpvaAOxdMpHqorFLjRGoBHm3lIXaURiUx7Q6mbEhERlGltm1QB28sfrYT5LzdvsFwaIzIDE3r1xLBTRyxYnRnFkFUb9aWMhZBdfBIgAuOzOqHj4YGw9rir+06svD4V/GIu5wjcjr6N5PvESosLERUVBRUKhVUKhXeeOMNTJgwodavZ48QNTYajQCplHf8kH6VKlR4Z+sZPNGhCQZ2aCJ2nEZDqdLg/e1nEXsxB/lldxdkDA90wRv9WqJHSxaZhsS7xv6kVquhUChga2uLsrIytG/fHqdOnYKbm1utXs9CiIxdXqkC6QXl6BTgInYUMmHL4lLw2a5LcLWT49A7fbhIp46KK6vwdexVfHv4BlSaux+7IX5O+N/ozvB24h5vhsChsT/JZDLY2t7dE0ahUEAQBO4iTCZDEATM+uksnlp2BBuO3xQ7DpmwV3oGYUiID1aN6cwiqA4crS3x3uBg/PhaV7jayQEAybeK0HVBLN7ekgyVSiNyQvMleiEUHx+PIUOGwMfHBxKJBNu3b69xzNKlSxEYGAhra2tERkbixIkTOr1HYWEhQkJC4Ofnh7feegvu7u56Sk8kri0Jt7D3QjZkUgkeYY8QGZClTIqvn+vE/87qqVNTFyS8H4U3+rWAhVQCQQB+OHULHT7ag+9PpIkdzyyJXgiVlZUhJCQES5cuvefzmzdvRnR0NObMmYPExESEhIRgwIAByMn5e8JZaGgo2rdvX+MrI+PuGg7Ozs5ITk5GamoqNm7ciOzs7Ab53ogM6dadcsz95QIAYHr/VmjbhEO31HCyiyux/xIn/taFRCLB9P6tkfTBY+jX5u48oYoqNWb9dBZjVp9A1j/uNiPDM6o5QhKJBNu2bcOwYcO0bZGRkQgPD0dMTAwAQKPRwN/fH1OnTsXMmTN1fo9Jkyahb9++GDly5D2fVygUUCgU2sfFxcXw9/fnHCEyKhqNgBdWHcfR6/kIa+qCH17rChknSFMDScsvx4hlh1FSqcL2yd1ZhNfTxcwiTNqQiNS8u/uV2VjKMLlPc7zSsxmsLWUip2u8TGKOkFKpREJCAqKiorRtUqkUUVFROHr0aK3OkZ2djZKSEgBAUVER4uPj0bp16/seP3/+fDg5OWm//P396/dNEBnA2qM3cPR6PmwsZVj0dAiLIGpQfi42aO/rhCB3O35Q60HbJk7YP6MPfprYDWFNXVBRpcZ/9lxBxw/3YP5vF6HRcP6QIRl1IZSXlwe1Wg0vL69q7V5eXsjKyqrVOW7evImePXsiJCQEPXv2xNSpU9GhQ4f7Hj9r1iwUFRVpv9LT0+v1PRDpW3pBORb8fgkA8O7gtgjkLuHUwKRSCb56thO2T+7OXer16JGmLtjyeld8NSoUdlYWUKo1+F/8dXSatxe7z9fuM490Z/JT/yMiIpCUlFTr462srGBlxYXoyHj5udhg9hPBOHQ1Dy9Gch8xEoeTrWW1xyWVVXCwtrzP0VRbEokET4b6oltzN0xcn4hTN++gqEKF175LQGsvB6x4KQxN3Vh86pNR9wi5u7tDJpPVmNycnZ0Nb29vkVIRiUsikeDFLk2x7MVHIJFwSIzE93PSbXRf8AcS0+6IHcVkeDhYY8vEbtj8ahd4Od794/xydgl6L4zDlI2JUKk5XKYvRl0IyeVyhIWFITY2Vtum0WgQGxuLrl27ipiMqOGVKlSorFJrH7MIImMgCAL2XMhGcaUK649xLSt9i2zmhuPvRuG9QW0gt5BCAPDrmUwMXnIIR67liR3PJIg+NFZaWopr165pH6empiIpKQmurq4ICAhAdHQ0xowZg86dOyMiIgKLFy9GWVkZxo0bJ2Jqoob3+a5L+ONSDuaP6ICeXJqfjIREIsGCER3Qyd8Z47oHiR3HZE14tDlGd2mKGVvOIP5qLi5nl+D5VcfRq5UHJvVujshmtdstgWoS/fb5uLg49OnTp0b7mDFjsGbNGgBATEwMFi5ciKysLISGhmLJkiWIjIxskHzcYoOMwaWsYgz66iA0ArBxQiS6NeeioETmqrBciS/3XsH642lQ/7ldR+emLlg++hG423O7jr9wrzE9YSFEYhOEu2sGHUnJx8D23lj2YpjYkYjuSxAErIi/jg6+TujWggW7IZ25VYjR35xAUUUVAEAqAZ6PCMDcJ9tBKjXqmS8NwiTWESIiYPf5LBxJyYfcQop3B7UVOw7RA609cgPzf7+EaZtOo7BcKXYck9bRzxmnZ0fh9V7NIJNKoBGA9cfT0OGjPfgp4ZbY8RoNFkJERqyySo2Pd14EALz2aDP4u9qKnIjowZ4ND0CInxOm928FJxveTm9oUqkUMwe2ReL7UejZ8m4PXJlCjegfk/HiquPILVE85AzEobGH4NAYiSnmj6v4z54r8Ha0xh8zenHXb2oU1BqBq52L5MytQkzekIj0OxUAAHsrC0zr1wJjuwVBbmFefR8cGiMyAVdzSgEAswa1YRFEjcY/iyClSoNLWcUipjEvHf2ccfCdvvju5Qh09HNCqUKFT3+7hNC5e7B47xWx4xkl9gg9BHuESGwJNwvwSIAL1w2iRievVIFX1p5Cal4Zdk7rAT8XDu02JI1GwJbEW5j7ywWUKlQAAFc7Ob54JgS9W3uKnM7w2CNEZCLCmrqyCKJGycHaAn/9rZ1eUCFyGvMjlUrwTGd/7J3+KDr6OQEACsqUGPvtSQz5+iBuF5aLnNA4sEfoIdgjRGL4Oek2ujZ3g6cD1wShxu3WnXIIAjjR3wgcuJyD6B+SkV92924+iQQY0ckXC0Z0hKUJzh9ijxBRI3UzvwzRPySj52f7kVVUKXYconrxc7GtVgTxb2/x9GrtiYTZ/TG9f0tYyiQQBGBr4m0M/+9hnLpRIHY80bAQIjIyy+JSoNYI6NrcDd5O7BEi03E+owhDYw7jRl6Z2FHM2hv9WuH07McQ1dYTtnIZzmUUY+Tyo5i8IRHJ6ea3cS4LISIjcruwAlsT7y6ENrVvC5HTEOnXgt8v4eztIny884LYUcyevbUFVo0JR/zbfTAq3B8SCbDzbCaeXHoEL6w8ZlaLYbIQIjIiy+NSUKUW0K25G8Kauoodh0ivFo4MwYhOvlg4MkTsKPQnd3srLHiqI36e1B0O1neX6Dicko+wj/fhk50XoNFoRE5oeCyEiIxEdnElNp9KBwBM7dtS5DRE+uftZI0vng2Fi51c7Cj0Lx39nZH8QX+M6RYIqeTuopgrD6YidO5e/HY2U+x4BsVCiMhI/O/AdShVGnRu6oIuzdgbRKbv0NU8XMspETsG/UkqleKjoe1w4r1+iAh0AQAUV6owaUMinl95DHfKTHO4jIUQkZGQW0hhZSHF1H4tuW4Qmbxtp29h9OrjmLQhEeVKldhx6B/c7a3xw+vd8MNrXeDlaAUAOJKSj97/icPaIzegUpvWcBnXEXoIriNEDSm/VAFXOzkLITJ5uSUKDFpyEP3aeOLDoe1gbSkTOxLdx69nMhDzxzVcyrrbe+dqJ8fEXs0w4dHmIid7sNp+frMQeggWQkREhpFXqoC7vZXYMagWVGoNvj+Zjs9+v6TdrsPLwQpLnuuEyGZuIqe7Ny6oSNRI/HomA6fTzG/tDqJ/F0HFlVUiJaGHsZBJMbpLU+yc1gNtvB0AANklCjy74hhGLjuCvNLGu/grCyEiEVWpNZj7ywUM/+8RHLqaJ3YcIlFUVqkxc+sZDIs5rO1tIOPU1M0Ou958FP97MQxONpYAgFM37yDik1i8v+1so7zdnoUQkYhiL+Ygp0QBNzs5IoJ4pxiZp3KlGgeu5CI1vwwHr+SKHYdqYUB7b5yeHYWJvZpDJpVAIwDrj6fhqWVHcfZWkdjxdMI5Qg/BOUJkSC+tPoH4K7mY2Ls53nm8jdhxiESTcPMOypUq9GzpIXYU0lFRuRJTvj+N49cLoFRrIJEAT4f5YULPZmjp5SBaLk6W1hMWQmQoafnleHThfkgkwIEZfRDgxt25iajxyiqqxILfL2J7Uoa2rV8bT3w1qhPs/1y1uiFxsjSRkdt4Ig0A8GhLDxZBRP9QWK7E+9vPooSTpxsVbydrLB7VCVsndoWD1d3CJ/ZSDjrN24OvYq+InO7+WAgRiUChUuPHP7fTeCEyQOQ0RMbl1e8SsP5YGmZvPyd2FKqDsKauOD27P54O84NEAlSpBXy59yrC5u3Fgcs5YsergYUQkQjS8sthZSGFt6M1+rbxFDsOkVGZNbANWnra45WezcSOQnVkYSHFwqdDcOidPujge3dYKr9MiTHfnsQLq44bVW8f5wg9BOcIkaGoNQLSCsoR5G4ndhQio6PWCJBJucK6qYi7nIPoH5JR8Od+Ze72crw9oA1GhvlBaqDrzMnSesJCiIhIXLfulMPB2lK7bg01XmuP3MDaIzdwPa8MAODrbIM3o1ri6c7+en8vTpYmMlIXM4tRZWKbFhIZSvyVXAz66iDe3pIM/t3e+I3pFohdbz6K9wa1ha1chtuFFXhryxnsF3HuEAshogZUoVTj2f8dRbcFf+BmfpnYcYiMnrOtJSqq1MgqVqC4kqtOmwK5hRQTHm2GbZO6oXNTF7TyskePFu6i5Wn4G/uJzNivZzJQXKmCn40l/F14yzzRw3T0c8aGV7og1N8Zcgv+7W5KWns7YsvEbqisUsNSJt61ZSFE1IA2HL+7dtDzkQEGmyBIZGr+vf2MIAiQSPjzYyqsLWWivj/La6IGciGjGEnphbCUSfB0mP4nBhKZOkEQsOH4TUzblMT5QqQ3LISIGsiO5LvLzvdr4wUPByuR0xA1PukFFfhoxwX8kpyBXeeyxI5DJoJDY0QNQBAE7Dx7txAaEuIjchqixinAzRazhwSjXKHCgHbeYschE8FCiKgBXMgsRnpBBawtpejThrtrE9XV6C5NxY5AJoaFEFEDCG7iiJ3TeuBaTils5fyxI9IHtUbAr2cyMDTEh5Onqc74G5moAUgkErTzcUI7HyexoxCZBEEQ8Mrak9h/ORe5JQruS0Z1xsnSRETU6EgkEvRt6wVbuYw3H1C9sEeIyMAW77uCtPxyjO0eiI5+zmLHITIZL0YGoF8bT/g424gdhRox9ggRGZAgCPjhZDp+On0bGYWVYschMikSiaRaEVRZpeb6QqQzFkJEBnQ6vRAZRZWwk8vQuzXvFiMylOu5pRi29DBWH74hdhRqZFgIERnQzjOZAIB+bb1EX0aeyJQdScnHpawSrIy/jgqlWuw41IhwjhCRgWg0An4/e7cQGtyxichpiEzbC5EBKKqowsgwP9jI+UcH1R4LISID+eewWK9WHBYjMiSJRILJfVqIHYMaIQ6NERnIX8Ni/YM5LEbU0BLT7uCHU+lix6BGgD1CRAYS6G6LVl72GNyRe4sRNaQLGcV4ZvlRSCRAG28HLltBD8RCiMhAXuoaiJe6BvJ2XqIG1raJAx5r5wWJRIIgdzux45CRYyFEZGDcA4moYUkkEnzxTCisLKT8+aOH4hwhIj3TaATsOpfFW3iJRGRtKatWBN0urBAxDRkzFkJEenYxqxivr09Az8/3Q6PhsBiRmFRqDT7ZeQF9FsYhOb1Q7DhkhFgIEelZ/JU8AECInxOkUnbLE4lJJpUgvaACSrUGh1PyxI5DRohzhIj07ODVXADAo1w7iEh0EokEn43siKc7+6FfWy+x45ARYo8QkR6VK1U4deMOAKBnS3eR0xARADjZWLIIovsyi0IoMDAQHTt2RGhoKPr06SN2HDJhx68XQKnWwM/FhrftEhmhUoUKM35MxrnbRWJHISNhNkNjR44cgb29vdgxyMQduHJ3WKxnSw/etktkhP6z+zK2JNzC6bQ72P3mo7CQmUV/AD2A2RRCRA3hr/lBvVpxWIzIGL0Z1RLnM4rwzuNtWAQRACMYGouPj8eQIUPg4+MDiUSC7du31zhm6dKlCAwMhLW1NSIjI3HixAmd3kMikaBXr14IDw/Hhg0b9JScqKYVL3XGnCHB6NqchRCRMXK2leOH17qic6Cr2FHISIjeI1RWVoaQkBCMHz8eI0aMqPH85s2bER0djeXLlyMyMhKLFy/GgAEDcPnyZXh6egIAQkNDoVKparx2z5498PHxwaFDh+Dr64vMzExERUWhQ4cO6Nixo8G/NzI/zT3s0dyDQ7BExuyfw9a5JQoUVVShhSd/bs2VRDCijZAkEgm2bduGYcOGadsiIyMRHh6OmJgYAIBGo4G/vz+mTp2KmTNn6vweb731Ftq1a4exY8fe83mFQgGFQqF9XFRUhICAAKSnp8PR0VHn9yMiIuN0/nYRJm9MhL2VBTa91hX2VqL3DZAeFRcXw9/fH4WFhXBycrrvcUZ91ZVKJRISEjBr1ixtm1QqRVRUFI4ePVqrc5SVlUGj0cDBwQGlpaX4448/8Mwzz9z3+Pnz5+Ojjz6q0e7v76/7N0BERI2C72yxE5ChlJSUNN5CKC8vD2q1Gl5e1dd/8PLywqVLl2p1juzsbAwfPhwAoFarMWHCBISHh9/3+FmzZiE6Olr7WKPRoKCgAG5ubtru1PDwcJw8efKer7/fc/9u/6tSNZaepgd9Tw15Pl1eV5tjH3ZMba/X/dp5Hev/Ol7Hmngda/fcvdqM6VryOtb+OUP8TAqCgJKSEvj4+DzwOKMuhPShWbNmSE5OrvXxVlZWsLKyqtbm7Oxc7bFMJrvvhbnfc/drd3R0FP2HFXjw99SQ59PldbU59mHH6Hq9eB31/zpex5p4HWv33IOON4ZryetY++cM9TP5oJ6gv4h+19iDuLu7QyaTITs7u1p7dnY2vL29RUoFTJ48WefnHvQaY6DvfHU9ny6vq82xDztG1+vF66j/1/E61sTrWLvneB31/zoxrmNt39dQGsVk6YiICHz99dcA7g5VBQQEYMqUKXWaLG0siouL4eTkhKKiItH/aqG643U0DbyOpoPX0jQ05HUUfWistLQU165d0z5OTU1FUlISXF1dERAQgOjoaIwZMwadO3dGREQEFi9ejLKyMowbN07E1PVnZWWFOXPm1BiGo8aF19E08DqaDl5L09CQ11H0HqG4uLh77v81ZswYrFmzBgAQExODhQsXIisrC6GhoViyZAkiIyMbOCkRERGZGtELISIiIiKxGPVkaSIiIiJDYiFEREREZouFEBEREZktFkJERERktlgINQLp6eno3bs3goOD0bFjR/z4449iR6I6Gj58OFxcXDBy5Eixo5AOfv31V7Ru3RotW7bEqlWrxI5DdcSfv8bPEJ+HvGusEcjMzER2djZCQ0ORlZWFsLAwXLlyBXZ2dmJHIx3FxcWhpKQEa9euxZYtW8SOQ7WgUqkQHByM/fv3w8nJCWFhYThy5Ajc3NzEjkY64s9f42eIz0P2CDUCTZo0QWhoKADA29sb7u7uKCgoEDcU1Unv3r3h4OAgdgzSwYkTJ9CuXTv4+vrC3t4eAwcOxJ49e8SORXXAn7/GzxCfhyyE9CA+Ph5DhgyBj48PJBIJtm/fXuOYpUuXIjAwENbW1oiMjMSJEyfq9F4JCQlQq9Xw9/evZ2r6t4a8jtRw6ntdMzIy4Ovrq33s6+uL27dvN0R0+gf+fJoGfV5HfX0eshDSg7KyMoSEhGDp0qX3fH7z5s2Ijo7GnDlzkJiYiJCQEAwYMAA5OTnaY0JDQ9G+ffsaXxkZGdpjCgoK8NJLL2HFihUG/57MUUNdR2pY+riuJD5eR9Ogr+uo189DgfQKgLBt27ZqbREREcLkyZO1j9VqteDj4yPMnz+/1uetrKwUevbsKaxbt05fUekBDHUdBUEQ9u/fLzz11FP6iEk6qst1PXz4sDBs2DDt82+88YawYcOGBslL91afn0/+/BmPul5HfX8eskfIwJRKJRISEhAVFaVtk0qliIqKwtGjR2t1DkEQMHbsWPTt2xejR482VFR6AH1cRzI+tbmuEREROHfuHG7fvo3S0lL8/vvvGDBggFiR6R7482kaanMdDfF5yELIwPLy8qBWq+Hl5VWt3cvLC1lZWbU6x+HDh7F582Zs374doaGhCA0NxdmzZw0Rl+5DH9cRAKKiovD000/jt99+g5+fH39Ji6w219XCwgKLFi1Cnz59EBoaiv/7v//jHWNGprY/n/z5M261uY6G+Dy0qNerqUH06NEDGo1G7BikB/v27RM7AtXB0KFDMXToULFjUD3x56/xM8TnIXuEDMzd3R0ymQzZ2dnV2rOzs+Ht7S1SKtIVr6Np4nU1DbyOpkGs68hCyMDkcjnCwsIQGxurbdNoNIiNjUXXrl1FTEa64HU0TbyupoHX0TSIdR05NKYHpaWluHbtmvZxamoqkpKS4OrqioCAAERHR2PMmDHo3LkzIiIisHjxYpSVlWHcuHEipqZ/43U0TbyupoHX0TQY5XXUy71nZm7//v0CgBpfY8aM0R7z9ddfCwEBAYJcLhciIiKEY8eOiReY7onX0TTxupoGXkfTYIzXkXuNERERkdniHCEiIiIyWyyEiIiIyGyxECIiIiKzxUKIiIiIzBYLISIiIjJbLISIiIjIbLEQIiIiIrPFQoiIiIjMFgshIiIiMlsshIiIiMhssRAiIqPXu3dvvPnmm2LHqDND5Z8xYwaGDRum9/MSmRMWQkRUa2PHjoVEIqnx9c/dpA3hp59+wrx58wxy7tzcXEycOBEBAQGwsrKCt7c3BgwYgMOHDxvk/fQpKSkJoaGhYscgatQsxA5ARI3L448/jm+//bZam4eHh8HeT6lUwtXVVS/nkcvlNdqfeuopKJVKrF27Fs2aNUN2djZiY2ORn59f7/c0tOTkZEyZMkXsGESNGnuEiEgnf/Wa/PNLJpNBoVBg2rRp8PT0hLW1NXr06IGTJ09We21gYCAWL15crS00NBQffvih9nHv3r0xZcoUvPnmm3B3d8eAAQNqDC1pNBrMnz8fQUFBsLGxQUhICLZs2VLtvPc6z78VFhbi4MGD+Oyzz9CnTx80bdoUERERmDVrFoYOHVrt/T7//HO0aNECVlZWCAgIwCeffAIA2LVrF3r06AFnZ2e4ubnhiSeeQEpKygP/DWuT/2Fu3bqFvLw8bY9QYWEhhgwZgh49eiArK0uncxGZMxZCRKQXb7/9NrZu3Yq1a9ciMTERLVq0wIABA1BQUKDzudauXQu5XI7Dhw9j+fLlNZ6fP38+1q1bh+XLl+P8+fOYPn06XnzxRRw4cECn89jb28Pe3h7bt2+HQqG4b55Zs2ZhwYIFmD17Ni5cuICNGzfCy8sLAFBWVobo6GicOnUKsbGxkEqlGD58ODQazX3PV9v8D5KUlARnZ2cEBgbi7NmzCA8Ph6+vL/bv3w9vb+9an4fI7AlERLU0ZswYQSaTCXZ2dtqvkSNHCqWlpYKlpaWwYcMG7bFKpVLw8fERPv/8c21b06ZNhS+//LLaOUNCQoQ5c+ZoH/fq1Uvo1KlTtWN69eolvPHGG4IgCEJlZaVga2srHDlypNoxL7/8svDcc8898Dz3smXLFsHFxUWwtrYWunXrJsyaNUtITk7WPl9cXCxYWVkJK1eufOi5BEEQcnNzBQDC2bNn65X/YebNmyf06tVL2LBhg+Di4iKsWLGi1q8lor9xjhAR6aRPnz5YtmyZ9rGdnR1SUlJQVVWF7t27a9stLS0RERGBixcv6vweYWFh933u2rVrKC8vR//+/au1K5VKdOrUqdbn+ctTTz2FwYMH4+DBgzh27Bh+//13fP7551i1ahXGjh2LixcvQqFQoF+/fvd8/dWrV/HBBx/g+PHjyMvL0/YEpaWloX379vXK/yBJSUk4c+YMpkyZgp07d6Jr1661fi0R/Y2FEBHpxM7ODi1atKjWlpubW6vXSqVSCIJQra2qquqe73E/paWlAICdO3fC19e32nNWVla1Ps8/WVtbo3///ujfvz9mz56NV155BXPmzMHYsWNhY2PzwNcOGTIETZs2xcqVK+Hj4wONRoP27dtDqVTWO/+DJCUlYcSIEdi4cSMKCwtr/Toiqo5zhIio3po3b66di/OXqqoqnDx5EsHBwdo2Dw8PZGZmah8XFxcjNTVVp/cKDg6GlZUV0tLS0KJFi2pf/v7+9f9m/nyPsrIyAEDLli1hY2OD2NjYGsfl5+fj8uXLeP/999GvXz+0bdsWd+7cMXj+kpISXL9+HZMnT0ZMTAxGjRqF8+fP6/6NEhF7hIio/uzs7DBx4kS89dZbcHV1RUBAAD7//HOUl5fj5Zdf1h7Xt29frFmzBkOGDIGzszM++OADyGQynd7LwcEBM2bMwPTp06HRaNCjRw8UFRXh8OHDcHR0xJgxY2p9rvz8fDz99NMYP348OnbsCAcHB5w6dQqff/45nnzySQB3e4veeecdvP3225DL5ejevTtyc3Nx/vx5jBs3Dm5ublixYgWaNGmCtLQ0zJw50+D5k5OTIZPJEBwcjE6dOuHcuXMYMmQITpw4AXd391p//0TEQoiI9GTBggXQaDQYPXo0SkpK0LlzZ+zevRsuLi7aY2bNmoXU1FQ88cQTcHJywrx583TuEQKAefPmwcPDA/Pnz8f169fh7OyMRx55BO+++65O57G3t0dkZCS+/PJL7Twnf39/TJgwodq5Zs+eDQsLC3zwwQfIyMhAkyZN8Prrr0MqlWLTpk2YNm0a2rdvj9atW2PJkiXo3bt3vfKvWbMG48aNqzGM+JekpCS0adNGO5S2cOFCXLx4ESNGjMC+ffvuuV4SEd2bRLjfTxoREYlizpw5OHDgAOLi4sSOQmTy2CNERGRkfv/9d8TExIgdg8gssEeIiIiIzBbvGiMiIiKzxUKIiIiIzBYLISIiIjJbLISIiIjIbLEQIiIiIrPFQoiIiIjMFgshIiIiMlsshIiIiMhssRAiIiIis8VCiIiIiMwWCyEiIiIyW/8P9SjR7XxH8TAAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(hm.k_hm, hm.power_1h_auto_tracer, ls=\"--\", color=\"C0\", label=\"1halo\")\n", "plt.plot(hm.k_hm, hm.power_2h_auto_tracer, ls=\":\", color=\"C0\", label=\"2halo\")\n", "plt.plot(hm.k_hm, hm.power_auto_tracer, color=\"C0\", label=\"full\")\n", "plt.legend()\n", "plt.xscale(\"log\")\n", "plt.yscale(\"log\")\n", "plt.ylim(\n", " 1e-5,\n", ")\n", "plt.xlabel(\"Fourier Scale, $k$\")\n", "plt.ylabel(\"Auto Power Spectrum\");" ] } ], "metadata": { "kernelspec": { "display_name": ".venv", "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.12.4" } }, "nbformat": 4, "nbformat_minor": 4 }