{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Customised extensions with `halomod`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this tutorial, we use the existing infrastructure of `halomod` and plug in a new type of tracer, namely HI using the model from [arxiv:2010.07985](https://arxiv.org/abs/2010.07985). This model requires three additions: a new density profile for HI; a new concentration-mass relation for HI; and finally a new HI HOD. \n", "\n", "`halomod` is extremely flexible. You can add custom models for *any* \"component\" (eg. mass functions, density profiles, bias functions, mass definitions,...). \n", "However, most likely you'll need to add something to build a new type of tracer, which is the case here." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's import a few basic things first:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import hmf\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "import halomod\n", "from halomod import TracerHaloModel" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using halomod v0.1.dev684+gc633d1a77 and hmf v3.5.2\n" ] } ], "source": [ "print(f\"Using halomod v{halomod.__version__} and hmf v{hmf.__version__}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Creating a new density profile" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The HI density profile used here is:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\rho_{\\rm HI} = \\rho_s \\bigg(\\frac{r_s}{r}\\bigg)^b {\\rm exp}\\bigg[-a \\frac{r}{r_s}\\bigg]\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that all the infrastructure has been set up by the `Profile` or `ProfileInf` class, depending on if you truncate the halos or not. And all you need to modify is the `_f` function, and optionally its integration `_h` and its Fourier transform `_p` (see documentation for [`profiles.py`](https://halomod.readthedocs.io/en/docs/_autosummary/halomod.profiles.html#module-halomod.profiles)). Note that these latter two are useful to speed up calculations, but are in principle optional and will be calculated numerically." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `_f` function is:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f(x) = \\frac{1}{x^b}{\\rm exp}\\big[-ax\\big]\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The integration in this case is:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "h = \\Gamma(3-b)\\times a^{b-3}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $\\Gamma$ is the Gamma function" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Fourier Transformed profile is:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p(K)= {\\rm tan}^{-1}(K/a)/K,\\,b=2\n", "$$\n", "$$\n", "\\begin{split}\n", "p(K)=&\\frac{-1}{1+(K/a)^2}\\Bigg(\\bigg(1+K^2/a^2\\bigg)^{b/2}\\times\\\\\n", "&\\Gamma(2-b)\\sin\\bigg[(b-2){\\rm arctan}\\big[K/a\\big]\\bigg]\\Bigg),\\,b\\ne2\n", "\\end{split}\n", "$$" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "from scipy.special import gamma\n", "\n", "from halomod.profiles import ProfileInf" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "class PowerLawWithExpCut(ProfileInf):\n", " \"\"\"\n", " A simple power law with exponential cut-off, assuming f(x)=1/x**b * exp[-ax].\n", " \"\"\"\n", "\n", " _defaults = ProfileInf._defaults | {\"a\": 0.049, \"b\": 2.248}\n", "\n", " def _f(self, x):\n", " return 1.0 / (x ** self.params[\"b\"]) * np.exp(-self.params[\"a\"] * x)\n", "\n", " def _h(self, c=None):\n", " return gamma(3 - self.params[\"b\"]) * self.params[\"a\"] ** (self.params[\"b\"] - 3)\n", "\n", " def _p(self, K, c=None):\n", " b = self.params[\"b\"]\n", " a = self.params[\"a\"]\n", " if b == 2:\n", " return np.arctan(K / a) / K\n", " else:\n", " return (\n", " -1\n", " / K\n", " * ((a**2 + K**2) ** (b / 2 - 1) * gamma(2 - b) * np.sin((b - 2) * np.arctan(K / a)))\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At bare minimum, you must specify `_f` which is just the density profile itself, and all the integration and Fourier transformation will be done numerically. However, that is very inefficient so you should always find analytical expression and specify if you can." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's plug it into a halo model:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "hm = TracerHaloModel(tracer_profile_model=PowerLawWithExpCut, transfer_model=\"EH\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that as of v2 of `halomod`, defining a new model for any particular component will automatically add it to a registry of 'plugins', and you may then construct the overall model using a string reference to the model (i.e. we could have passed `tracer_profile_model=\"PowerLawWithExpCut\"`). " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And see the profile for a halo of mass $10^{10} M_\\odot h^{-1}$ in Fourier space:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGhCAYAAACHw3XjAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAANZZJREFUeJzt3Xl8VOWh//HvmclGVgghCQlhkzUBEgsBF1DQKKKCS1162ypiy22tWr1RW5eWRa3carW2Ovf6U2vFVm9RW5FFqYJYKosoGJDdyBYICYSYnWwz8/sjGjcCCZnJM3Pm83698oKZOTPzzfMaZr7Mec55LK/X6xUAAIBNOUwHAAAA8CfKDgAAsDXKDgAAsDXKDgAAsDXKDgAAsDXKDgAAsDXKDgAAsLUw0wF8zePxqLi4WHFxcbIsy3QcAADQDl6vV9XV1UpLS5PD4dvvYmxXdoqLi5WRkWE6BgAAOAVFRUXq06ePTx/TdmUnLi5OUstgxcfHG04DAADao6qqShkZGa2f475ku7Lzxa6r+Ph4yg4AAEHGH1NQmKAMAABsjbIDAABsjbIDAABsjbIDAABsjbIDAABsjbIDAABsjbIDAABsjbIDAABsjbIDAABsLSDLzpIlSzR06FANHjxYzz77rOk4AAAgiAXcchHNzc3Kz8/XypUrlZCQoNGjR+uKK65Qz549TUcDAABBKOC+2Vm/fr2ysrKUnp6u2NhYTZkyRW+99ZbpWAAAIEj5/JudVatW6ZFHHtGGDRt06NAhvfbaa7r88su/to3L5dIjjzyikpISZWdn64knntDYsWMlScXFxUpPT2/dNj09XQcPHuxwjv3ltYprch73NkttLzLWmfXHTnTfky1sdqJbT/i4J7znye57Eqf4vCcbwxP/rid43BM/bKfG6YS/6wluc1iWwhyWHJYlp8OSw/LPInYAgFPn87JTW1ur7Oxs3Xjjjbryyiu/dfuCBQuUn5+vp556SuPGjdPjjz+uyZMna+fOnUpOTu7w8zU0NKihoaH1clVVlSTp4j+8J0dk9Kn/IsApcjosOS1LDocU5nDIYX1+ncMhp0NyWpaczi+2sRThdCgy3KlIp0OR4Q5FhjkUEeZQZJjzK39vuRwTGabYqDDFRYa1/D0yTHFRX/7ZPTpCTgdlCwC+yudlZ8qUKZoyZUqbtz/22GOaOXOmZsyYIUl66qmntHTpUj333HO6++67lZaW9rVvcg4ePNj6rc/xzJs3T3Pnzv3W9dERDjkjvv3Njvck+b0n2cB7gkc4+X1P9uTmntt7kgc40a0ne+5Q4/Z45ZZXckuSp0uf27Kk7t3ClRgToZ4xkS1/xkYorXs3pXWPUnr3aKV1j1JKfJTCnQG3FxsA/MLynuxTrjMPbllf243V2Nio6Ohovfrqq1/btTV9+nRVVFTo9ddfV3Nzs4YPH6533323dYLymjVr2pygfLxvdjIyMlRZWan4+Hh//WrwoZMWrRPc7M8Sd/Ln9srrbSk3zR6vPB6v3F5vS9n56o+35bbmzy97vrJNs8erJrdHjc0eNTR71NDs/vLvTV+/XN/kVk2DW7UNzappaFZ1Q3PL3+tbLtc0NJ/kt/mSw5JS46M0oFeMBvWK1aDkWJ2WHKuhKXHqGRvZ7scBAF+pqqpSQkKCXz6/u/RorLKyMrndbqWkpHzt+pSUFO3YsaMlUFiYHn30UU2aNEkej0e/+MUvTngkVmRkpCIjeXMOZiedz9SpvTKhs0un2e3RZ3VNKq9t1NHaBpXXNqq8tlFl1Q06WFGv4opjKq48pkMV9Wp0e1RcWa/iynqtLjz6tcdJ795No/okaFSf7srOSFBORndFRwTcgZsA0G4B+Q42bdo0TZs2zXQMIKiEOR3qFRepXnGRkuLa3M7j8aqstkFF5cf06ZEafXq4RoWHa1R4pEb7y+t0sOKYDlYc05tbSiRJ4U5Lp/ftobNPS9L4wT01qk93doEBCCpdWnaSkpLkdDpVWlr6tetLS0uVmpralVGAkOVwWEqOi1JyXJRG9+vxtduq65v08cFKbT5Qqc0HKlSwv0LFlfVav6dc6/eU6/fLpbjIMJ0/PFkXj+ytc4b0UlT48Y96BIBA0aVlJyIiQqNHj9aKFSta5+x4PB6tWLFCt9xyS1dGAXAccVHhOuu0JJ11WpKkljlP+8vr9F5hmVYXlmnNp0dVUdekhQXFWlhQrNjIMOUNT9bU7DSdO6SXwvjGB0AA8nnZqampUWFhYevlPXv2qKCgQImJierbt6/y8/M1ffp0jRkzRmPHjtXjjz+u2tra1qOzAAQOy7LUr2eM+vWM0Q/G9ZPH49VHRZ9p6eYSvbnlkA5V1rcWn9T4KF2Tm6FrczOU3r2b6egA0MrnR2O9++67mjRp0reunz59up5//nlJ0pNPPtl6UsGcnBz98Y9/1Lhx4zr1vC6XSy6XS263W7t27eJoLMDPWopPhZZsLtbCjw7qs7omSS3nFLp4ZG/95JyBGpGeYDglgGDhz6Ox/HrouQn+HCwAx9fQ7NZbW0v14vv7tG53eev1Zw/qqf/KG6Ix/RMNpgMQDCg7HUDZAczacrBSz/57txZvPiS3p+XtZeLQXrrjgqEa2YdvegAcH2WnAyg7QGA48FmdXCsL9fKHB1pLz5Wnp+uXU4YpJT7KcDoAgYay0wGUHSCw7C2r1R9WfKLXPmpZBiY6wqlbzhukH48fqIgwjt4C0IKy0wGUHSAwbSqq0JzFW/XR/gpJ0tCUOP32qlHKyehuNBeAwODPz2/+WwWgS2RndNfff3qWHrsmWz1jIrSztFpX/s9qPbBkm441uk3HA2BjlB0AXcbhsHTld/ro7fxzdcXp6fJ4pT+9t0dTn3xP24qrTMcDYFO2KTsul0uZmZnKzc01HQXASSTGROj31+bozzfkqldcpAoP1+hy12o9++/d8nhstWcdQABgzg4Ao47WNOiXf/9Yy7e3rJmXNzxFj16TrYRu4YaTAehKzNkBYFs9YyP1zPWj9cDlIxQR5tDy7aW67Mn3tLOk2nQ0ADZB2QFgnGVZuu6Mfnr1p2cqvXs37T1ap8tdq7V08yHT0QDYAGUHQMAY1ae7Ft86XhMGJ+lYk1s3v7RRrpWFstnedgBdjLIDIKAkxkTo+RljdePZAyRJj/xzp+58ZbMamz2GkwEIVpQdAAHH6bA0a2qmHrh8hJwOS3/feEAznl+v2oZm09EABCHKDoCAdd0Z/fTcDbmKiXBqdeFR/eDZ91VR12g6FoAgQ9kBENDOHdJLL848Q92jw1VQVKFr/t9alVbVm44FIIhQdgAEvJyM7nr5J2cqJT5Su0prdNVTa7TvaK3pWACChG3KDmdQBuxtSEqcXv3pWerXM1pF5cd01VNrVXi4xnQsAEGAMygDCCqHq+t1/Z/Wa0dJtZLjIrXgJ2dqQFKM6VgAOokzKAPA55LjovTSzDM0NCVOh6sb9P1n1qmovM50LAABjLIDIOgkxkToxZnjdFqvGB2qrNf3nl6ngxXHTMcCEKAoOwCCUlJspP5v5hkakBSjgxXH9B9Pr1NJJUdpAfg2yg6AoJUcH6WXZo5T38Ro7S+v0w+eXafPajkPD4Cvo+wACGq9E7rppZnjlJYQpU+P1OrG+R+orpEzLQP4EmUHQNDr0yNa828cq4Ru4fpof4VueekjNbtZSwtAC8oOAFsYnBKn524Yo8gwh97ZcVj3vvYxq6UDkETZAWAjo/sl6snvf0cOS3r5wwN69K1dpiMBCAC2KTucQRmAJF2QmaKHrhgpSXpyZaH+um6f4UQATOMMygBs6Q/LP9Hvl++S02Hp+Rm5mjC4l+lIAE6AMygDQAf9/PxBuuL0dLk9Xv3sxY2sowWEMMoOAFuyLEv//d2RGt2vh6rrm/Wj+R9wDh4gRFF2ANhWZJhT/++60erTo5v2Ha3TT/66QY3NHJIOhBrKDgBbS4qN1HM35Co2Mkzr95TrPg5JB0IOZQeA7Q1JidOT3z9dDkt6ZcMBzV+z13QkAF2IsgMgJEwcmqx7Lx4uSXpg6Xa9v/uo4UQAugplB0DI+NH4AZqWnSa3x6ubX9qoQ5XHTEcC0AUoOwBChmVZ+u13R2lYapzKahr1079uVEOz23QsAH5G2QEQUrpFOPX0dWOU0C1cm4oqNGvhViYsAzZH2QEQcvr2jNYT/9EyYXnBh0V6af1+05EA+BFlB0BIOmdIL905eagkae6ibdpysNJwIgD+Ypuyw0KgADrqpnNPU97wZDW6Pbr5pY2qqm8yHQmAH7AQKICQVlHXqEv++J4OVhzTxSNT5fr+d2RZlulYQMhhIVAA8JPu0RF68vunK9xp6Y2PS/TC2n2mIwHwMcoOgJB3et8euntKywkHH1y6TZsPVJgNBMCnKDsAIOnGs/trclaKmtwtJxysPMb8HcAuKDsAoJYTDj58VbYyErupqPyYfvHqJs6/A9gEZQcAPpfQLVyu739HEU6H/rm1VH9evdd0JAA+QNkBgK8Y1ae77rukZf7OvDe3c/4dwAYoOwDwDdef2U8XZrbM37ntbx/pWCPrZwHBjLIDAN/wxYKhKfGR+vRIrR5Yus10JACdQNkBgOPoEROhx67JkWVJL72/X8u2lJiOBOAUUXYAoA1nD0rSf54zUJJ09z82q6Sy3nAiAKeCsgMAJ3DHBUM1Mj1BFXVNyn+5QB4Ph6MDwYayAwAnEBHm0B++l6Nu4U6t+fSonv73btORAHQQZQcATmJgr1jNnZYlSfrdP3eynAQQZCg7ANAOV4/po4tHpqrZ49VtfytQXWOz6UgA2sk2ZcflcikzM1O5ubmmowCwIcuyNO+KUeqdEKU9ZbWa98YO05EAtJPltdniL1VVVUpISFBlZaXi4+NNxwFgM6sLy/SDZ9+XJM2/cazOHdLLcCLAHvz5+W2bb3YAoCucPShJN5zVX5L0i1c3qaKu0WwgACdF2QGADvrlRcM0sFeMSqsaNOv1rabjADgJyg4AdFC3CKd+f02OnA5LizYVa/GmYtORAJwAZQcATkF2RnfdMmmQJOlXC7eotIqzKwOBirIDAKfolvMGaVSfBFUea9Jdr26WzY73AGyDsgMApyjc6dBj12QrMsyhVbuO6MX395uOBOA4KDsA0AmDkuP0y4uGSZJ+s3S79h2tNZwIwDdRdgCgk244q7/OHNhTx5rcuuvVzSwWCgQYyg4AdJLDYenhq0YpOsKp9XvK9Zd1+0xHAvAVlB0A8IGMxGjdM6Vld9Z/v7mD3VlAAKHsAICP/GBcv9bdWb9gdxYQMCg7AOAjX92d9T67s4CAQdkBAB/KSIzW3V/ZnbX/aJ3hRAAoOwDgYz8c109nDExs2Z31903szgIMo+wAgI85HJYe/m62uoU7tW53uf76PruzAJMoOwDgB317sjsLCBSUHQDwk+vO6KdxAxJV18juLMAkyg4A+InDYemRq77cnfUiu7MAI2xTdlwulzIzM5Wbm2s6CgC0+ururHlv7tCBz9idBXQ1y+v12up71aqqKiUkJKiyslLx8fGm4wCAPB6vrn16rT7Y+5nOHdJLz8/IlWVZpmMBAcWfn9+2+WYHAAKVw2Fp3pWjFOF06F+7jmhhwUHTkYCQQtkBgC4wKDlWt+UNliTdv3ibjtY0GE4EhA7KDgB0kf88Z6CG947XZ3VNmrt4m+k4QMig7ABAFwl3OvTb746Uw5IWbSrWOztKTUcCQgJlBwC60Kg+3fXjCQMlSfe9tkXV9U2GEwH2R9kBgC72X3lD1K9ntA5V1uvhZTtNxwFsj7IDAF2sW4RT864YKUn6y7p9+mBvueFEgL1RdgDAgLMGJenaMRmSpF/+fbPqm9yGEwH2RdkBAEPuvXi4esVFaveRWj35TqHpOIBtUXYAwJCE6HA9cFmWJOmpf32q7YeqDCcC7ImyAwAGXTSity7KSlWzx6tf/n2zmt0e05EA26HsAIBh91+WpfioMG0+UKk/r95rOg5gO5QdADAsOT5K910yXJL06Ns7tf8oK6MDvkTZAYAAcM2YDJ11Wk/VN3n0q9e3yOv1mo4E2AZlBwACgGVZevDyEYoIc2jVriNatKnYdCTANig7ABAgBvaK1S2TBkmSHliyTZV1LCUB+AJlBwACyE/PPU2DkmNVVtOo/1623XQcwBYoOwAQQCLCHHro86Uk/m99kdbvYSkJoLMoOwAQYMYOSNT3cluWkrj3tY/V2My5d4DOoOwAQAC6Z8pwJcVGqPBwjf7fvz41HQcIapQdAAhACdHh+vWlmZKkJ1YWaveRGsOJgOBF2QGAADUtO00TBiepsdmjXy3k3DvAqaLsAECAsixLv7l8pKLCHVrz6VH9Y+NB05GAoETZAYAA1rdntG47f4gk6cGl21Re22g4ERB8KDsAEOB+PGGAhqXG6bO6Jj30BufeATqKsgMAAS7c6dBDV46UZUmvbjigNZ+WmY4EBBXblB2Xy6XMzEzl5uaajgIAPvedvj30w3H9JEn3vbZF9U1uw4mA4GF5bTa9v6qqSgkJCaqsrFR8fLzpOADgM1X1Tcp79F86XN2gn58/WPkXDDEdCfAZf35+2+abHQCwu/iocM2ZliVJ+t93C1V4uNpwIiA4UHYAIIhMGZGq84clq8nt1b3/4Nw7QHtQdgAgiFiWpbmXZalbuFPr95br1Q0HTEcCAh5lBwCCTJ8e0bo9b7Akad6bO/QZ594BToiyAwBB6MbxAzQ0JU7ltY367bIdpuMAAY2yAwBBKNzp0INXjJAk/e2DIm3YV244ERC4KDsAEKRy+yfqmjF9JLWce6fJ7TGcCAhMlB0ACGJ3Txmu7tHh2lFSredX7zUdBwhIlB0ACGKJMRG6Z8owSdLvl+9SccUxw4mAwEPZAYAgd/XoDI3p10N1jW7dv3ib6ThAwKHsAECQczgsPXjFCDkdlpZtLdE7O0pNRwICCmUHAGxgWGq8fjR+gCRp1utbdayRhUKBL1B2AMAmbjt/sNISonTgs2N6cuUnpuMAAYOyAwA2ERMZptmfLxT69KrdLBQKfI6yAwA2cmFmSutCob9ayEKhgETZAQBbsSxLc6ZlKSrcoXW7y/XaRwdNRwKMo+wAgM1kJEbr5+e3LBT6m6XbVVHHQqEIbZQdALChH48fqMHJsTpa26iH/7nTdBzAKMoOANhQRJhDD1zeslDo/63fr437PzOcCDCHsgMANnXGwJ668jvp8npbFgptZqFQhCjKDgDY2L0XD1dCt3BtP1Sl+Wv3mY4DGEHZAQAbS4qN1C8valko9LG3dqqkst5wIqDrUXYAwOa+l5uh0/t2V22jWw8uZaFQhB7KDgDYnMNh6YHLRshhSUs2H9LqwjLTkYAuRdkBgBAwIj1B153RT5L069e3qLGZycoIHZQdAAgR+RcOVVJshHYfqdWz7+02HQfoMpQdAAgRCd3Cdc+U4ZKkJ1YU6mDFMcOJgK5B2QGAEHLld9KV27+HjjW59cBiJisjNFB2ACCEWJalBy4fIafD0rKtJXp352HTkQC/o+wAQIgZlhqvG87qL0mavWir6pvcZgMBfkbZAYAQdHveYCXHRWrf0To9vYrJyrA3yg4AhKC4qHDdd0nLZGXXykIVldcZTgT4D2UHAELUtOw0nTmwpxqaPZq7eKvpOIDfUHYAIERZlqX7L8tSmMPS8u2HtXxbqelIgF9QdgAghA1OidOPJgyQJM1ZzGRl2BNlBwBC3M/PG6zeCVE68Nkx/c/KQtNxAJ+j7ABAiIuJDNOvL82UJD31r93aW1ZrOBHgW5QdAICmjEjVhMFJanR7NHvRVnm9XtORAJ+h7AAAPp+sPEIRTof+teuI/rm1xHQkwGcoOwAASdKApBj95zkDJUn3L96musZmw4kA36DsAABa3TxpkNK7d1NxZb2eeIfJyrCHgCw7V1xxhXr06KGrrrrKdBQACCndIpyaPbVlsvKz/96twsM1hhMBnReQZee2227TCy+8YDoGAISkCzJTdN6wZDW5vZq9aAuTlRH0ArLsTJw4UXFxcaZjAEBIsixLc6ZmKSLModWFR7Vk8yHTkYBO6XDZWbVqlaZOnaq0tDRZlqWFCxd+axuXy6X+/fsrKipK48aN0/r1632RFQDQRfr2jNbPJp4mSXpw6TbVNDBZGcGrw2WntrZW2dnZcrlcx719wYIFys/P1+zZs7Vx40ZlZ2dr8uTJOnz4cOs2OTk5GjFixLd+iouLO/wLNDQ0qKqq6ms/AIDO++m5p6lfz2iVVjXoD8t3mY4DnLKwjt5hypQpmjJlSpu3P/bYY5o5c6ZmzJghSXrqqae0dOlSPffcc7r77rslSQUFBaeW9jjmzZunuXPn+uzxAAAtosKdmjMtSzP+/IGeW71XV43O0NBUphgg+Ph0zk5jY6M2bNigvLy8L5/A4VBeXp7Wrl3ry6dqdc8996iysrL1p6ioyC/PAwChaNLQZF2YmSK3x6tfv85kZQQnn5adsrIyud1upaSkfO36lJQUlZS0/2yceXl5uvrqq/XGG2+oT58+JyxKkZGRio+P/9oPAMB3Zk3NVFS4Q+v3lGthwUHTcYAOC8ijsZYvX64jR46orq5OBw4c0Jlnnmk6EgCErD49onXreYMlSb9ZukNV9U2GEwEd49Oyk5SUJKfTqdLS0q9dX1paqtTUVF8+FQCgC/14wgANTIpRWU2DHn/7E9NxgA7xadmJiIjQ6NGjtWLFitbrPB6PVqxYwbczABDEIsNaJitL0vy1e7WjhCNfETw6XHZqampUUFDQekTVnj17VFBQoP3790uS8vPz9cwzz2j+/Pnavn27brrpJtXW1rYenQUACE7nDOmli7JS5fZ4Nfv1rUxWRtDo8KHnH374oSZNmtR6OT8/X5I0ffp0Pf/887r22mt15MgRzZo1SyUlJcrJydGyZcu+NWkZABB8fnXpcL2767De31OuRZuKdVlOuulIwElZXptUc5fLJZfLJbfbrV27dqmyspIjswDAD55Y8YkefXuXUuIjteKOiYqN7PD/m4FvqaqqUkJCgl8+v21Tdr7gz8ECAEj1TW5NfnyV9h2t00/OGah7Lh5uOhJswJ+f3wF56DkAIHBFhTs1e2qmJOlP7+1R4eFqw4mAE6PsAAA67LxhKcobnqxmj1dzFm1jsjICGmUHAHBKZl2apYgwh94rLNObW9p/lnygq1F2AACnpG/PaP303NMkSQ8u2aa6xmbDiYDjo+wAAE7Zzyaepj49uqm4sl6ulYWm4wDHRdkBAJyyqHCnfn1py2TlZ1bt0Z6yWsOJgG+zTdlxuVzKzMxUbm6u6SgAEFIuzEzRuUN6qdHt0ZxFnFkZgYfz7AAAOm1PWa0m/36VGt0ePX3daF2YxeLP6BjOswMACGgDkmL04wkDJEn3L9mm+ia34UTAlyg7AACfuOW8QUpLiNKBz47pf9/91HQcoBVlBwDgE9ERYfrV55OV//dfn2r/0TrDiYAWlB0AgM9MGZGqswf1VGOzR/cv2Wo6DiCJsgMA8CHLsjR3WpbCHJaWbz+sd3aUmo4EUHYAAL41KDlOPxrfMll57mImK8M8yg4AwOduPX+wUuIjte9onZ79927TcRDibFN2OKkgAASO2Mgw3XvxcEnSkysLdeAzJivDHE4qCADwC6/Xq+89vU7v7ynXRVmpeuq60aYjIYBxUkEAQNCxLEv3XzZCToelZVtLtGrXEdOREKIoOwAAvxmaGqfpZ/aXJM1ZvFWNzR6zgRCSKDsAAL+6/YLBSoqN1O4jtfrTe3tMx0EIouwAAPwqPipc90wZJkl64p1PdKjymOFECDWUHQCA3135nXSN6ddDdY1u/WbpdtNxEGIoOwAAv7MsS3Mvy5LDkpZsPqQ1n5aZjoQQQtkBAHSJrLQE/fCMfpKk2a9vVZObycroGpQdAECXueOCoeoZE6FPDtdo/pq9puMgRFB2AABdJiE6XL+8qGWy8uPLP9HhqnrDiRAKbFN2WC4CAILDVaP7KCeju2oamjXvzR2m4yAEsFwEAKDLbT5Qoctcq+X1Si//5EyNHZBoOhIMY7kIAICtjOrTXd/L7StJmvX6FjUzWRl+RNkBABjxi8lD1T06XDtKqvXXdftMx4GNUXYAAEb0iInQXZOHSpIefXuXymoaDCeCXVF2AADGfC+3r0akx6u6vlmPLNtpOg5sirIDADDG6bA0d1qWJOnlDUUqKKowGwi2RNkBABg1ul+irvxOurxeafbrW+Tx2OogYQQAyg4AwLi7pwxTbGSYNh2o1KsbDpiOA5uh7AAAjEuOi9LteYMlSb9dtkOVx5oMJ4KdUHYAAAFh+ln9NSg5VkdrG/X48l2m48BGKDsAgIAQ7nRoztSWycovrN2nnSXVhhPBLig7AICAMX5wkqaMSJXb49XsRVtksxWNYIhtyg4LgQKAPdx3yXBFhTu0bne5ln58yHQc2AALgQIAAs4fln+i3y/fpd4JUVpxx7mKjggzHQl+xkKgAICQ8pNzByojsZsOVdbLtbLQdBwEOcoOACDgRIU79etLMiVJz6zao71ltYYTIZhRdgAAAemCzBSdM6SXGt0e3b9km+k4CGKUHQBAQLIsS7OnZircaemdHYf1zo5S05EQpCg7AICAdVqvWN04foAkae7ibapvchtOhGBE2QEABLRbzxus5LhI7Ttapz+9t8d0HAQhyg4AIKDFRobpvkuGS5KefKdQxRXHDCdCsKHsAAAC3rTsNI3tn6hjTW499MZ203EQZCg7AICAZ1mW5kzLksOSlmw+pDWflpmOhCBC2QEABIXMtHj98Ix+kqS5i7ap2e0xnAjBgrIDAAga+RcMUY/ocO0srdZf1u0zHQdBgrIDAAga3aMjdNfkYZKkx97epbKaBsOJEAwoOwCAoHJtboZGpMerur5ZDy/bYToOgoBtyo7L5VJmZqZyc3NNRwEA+JHTYWnutBGSpJc/PKCCogqzgRDwLK/X6zUdwpf8uUQ8ACBw3PHyJv194wFl90nQaz87Ww6HZToSOsGfn9+2+WYHABBafjllqGIjw7TpQKVe2VBkOg4CGGUHABCUkuOidHveYEnSw8t2qrKuyXAiBCrKDgAgaE0/q78GJ8fqaG2jfr98l+k4CFCUHQBA0Ap3OjRnWpYk6S/r9mlHSZXhRAhElB0AQFA7e1CSLh6ZKrfHq9mvb5XNjruBD1B2AABB775LMhUV7tD7e8q1ePMh03EQYCg7AICgl969m342cZAk6aGl21Xb0Gw4EQIJZQcAYAv/ec5A9U2MVklVvVwrC03HQQCh7AAAbCEq3KlfX5opSXrm37u1p6zWcCIECsoOAMA28oYn69whvdTk9ur+xVtNx0GAoOwAAGzDsizNnpqpcKellTuPaMX2UtOREAAoOwAAWxnYK1Y/Gj9QkjR38TbVN7kNJ4JplB0AgO3cet4gpcRHan95nZ79927TcWAYZQcAYDsxkWG69+LhkqQnVxbqYMUxw4lgEmUHAGBL07LTNHZAouqbPHpo6XbTcWAQZQcAYEuWZWnO1Cw5LGnpx4e0prDMdCQYQtkBANhWZlq8rjujnyRpzuKtanJ7DCeCCZQdAICt5V8wVIkxEdpVWqMX1u4zHQcG2KbsuFwuZWZmKjc313QUAEAASYgO112Th0qSHn97l45UNxhOhK5meb1er+kQvlRVVaWEhARVVlYqPj7edBwAQABwe7y64n9Wa/OBSl09uo8euTrbdCR8gz8/v23zzQ4AAG1xOizNnZYlSXplwwFt3P+Z4UToSpQdAEBIOL1vD101uo8kac6irfJ4bLVjAydA2QEAhIxfXjRMcZFh2nygUi9/WGQ6DroIZQcAEDJ6xUXq9guGSJIe/udOVdY1GU6ErkDZAQCElOvP7KfBybEqr23UY2/vNB0HXYCyAwAIKeFOR+tk5b+s26fth6oMJ4K/UXYAACHnrEFJumRkb3m80uzXt8pmZ2HBN1B2AAAh6d5Lhisq3KH1e8u1aFOx6TjwI8oOACAkpXfvppsnDpIkPfTGdtU2NBtOBH+h7AAAQtbMcwaqb2K0Sqsa9MQ7habjwE8oOwCAkBUV7tSsSzMlSX96b7d2H6kxnAj+QNkBAIS084cna9LQXmpye3X/km1MVrYhyg4AIKRZlqVfX5qpcKeld3ce0Yrth01Hgo9RdgAAIW9gr1j9aPxASdL9S7apvsltOBF8ibIDAICkW88bpJT4SO0vr9Of3ttjOg58iLIDAICkmMgw3TNluCTpyXcKVVxxzHAi+AplBwCAz12Wk6bc/j10rMmteW/uMB0HPkLZAQDgc5Zlac60LDksafGmYq3bfdR0JPgAZQcAgK/ISkvQf4ztK0mas2irmt0ew4nQWZQdAAC+4c4Lh6p7dLh2lFTrpfX7TcdBJ1F2AAD4hh4xEbrjwqGSpEff2qXy2kbDidAZlB0AAI7j+2P7anjveFUea9Ij/9xpOg46gbIDAMBxOB2W5k7LkiT97YP92nKw0nAinCrKDgAAbRg7IFGX5aTJ65VmL9rKullBirIDAMAJ3DNluKIjnNqw7zO99tFB03FwCig7AACcQGpClG45b5Akad6bO1Rd32Q4ETrKNmXH5XIpMzNTubm5pqMAAGzmR+MHaEBSjI5UN+jJdwpNx0EHWV6b7YCsqqpSQkKCKisrFR8fbzoOAMAmVu44rBnPf6Bwp6Vlt5+j03rFmo5kK/78/LbNNzsAAPjTpGHJOm9YsprcXs1dvI3JykGEsgMAQDvNujRTEU6HVu06ouXbD5uOg3ai7AAA0E79k2L04wkDJEkPLNmm+ia34URoD8oOAAAdcPOkQUqNj9L+8jo9s2q36ThoB8oOAAAdEBMZpnsuHiZJcr1bqOKKY4YT4WQoOwAAdNC07DSN7Z+o+iaPfvPGdtNxcBKUHQAAOsiyLM2ZliWHJS3dfEhrPi0zHQknQNkBAOAUZKbF6wfj+kmS5i7apma3x3AitIWyAwDAKbrjwiHqHh2unaXV+uu6fabjoA2UHQAATlH36AjdeeFQSdJjb+/S0ZoGw4lwPJQdAAA64T/G9lVWWryq6pv1u7d2mo6D46DsAADQCU6HpbnTsiRJf/ugSJsPVJgNhG+h7AAA0Elj+ifq8pw0eb3S7EVb5fGwblYgoewAAOAD91w8XDERTn20v0L/+Oig6Tj4CsoOAAA+kBIfpVvPHyxJ+u83d6i6vslwInyBsgMAgI/MOLu/BiTFqKymQX9c8YnpOPgcZQcAAB+JDHNq1tRMSdKfV+9V4eFqw4kgUXYAAPCpSUOTlTc8Wc0er+Yu3iavl8nKplF2AADwsV9dkqkIp0P//qRMb20rNR0n5FF2AADwsf5JMZp5zgBJ0gNLtqm+yW04UWij7AAA4Ac3Txqk1PgoHfjsmJ5etdt0nJBG2QEAwA+iI8J07yXDJUn/826hDnxWZzhR6KLsAADgJ1NH9da4AYmqb/LooTe2m44Tsig7AAD4iWVZmjMtSw5LeuPjEq0pLDMdKSRRdgAA8KPhveP1wzP6SZLmLN6qJrfHcKLQQ9kBAMDP8i8Yoh7R4dpVWqO/rN1nOk7IoewAAOBn3aMjdNfkYZKk3y/fpbKaBsOJQgtlBwCALnBtboZGpMerur5ZjyzbaTpOSKHsAADQBZwOS3OnZUmSXt5QpE1FFWYDhRDKDgAAXWR0v0RdeXq6vF5p9qKt8nhYN6srUHYAAOhCd08ZppgIpwqKKvT3jQdMxwkJlB0AALpQcnyUfn7+YEnSb5ftUFV9k+FE9kfZAQCgi804e4AGJsWorKZRf1j+iek4tkfZAQCgi0WEOTRraqYkaf6avfqktNpwInuj7AAAYMDEocnKG56iZo9Xcxdvk9fLZGV/oewAAGDIrEszFRHm0HuFZfrn1hLTcWyLsgMAgCF9e0brPycMlCQ9sGS76pvchhPZE2UHAACDfjbpNKUlROlgxTE99a9PTcexJcoOAAAGRUeE6d5LhkuS/vfdT1VUXmc4kf1QdgAAMOySkb11xsBENTR79NAb203HsR3KDgAAhlmWpTnTsuR0WHpzS4lWF5aZjmQrAVd2ioqKNHHiRGVmZmrUqFF65ZVXTEcCAMDvhqXG67oz+klqWTerye0xnMg+Aq7shIWF6fHHH9e2bdv01ltv6fbbb1dtba3pWAAA+N1/5Q1RYkyECg/X6IW1+0zHsY2AKzu9e/dWTk6OJCk1NVVJSUkqLy83GwoAgC6QEB2uuyYPlSQ9/vYuHaluMJzIHjpcdlatWqWpU6cqLS1NlmVp4cKF39rG5XKpf//+ioqK0rhx47R+/fpTCrdhwwa53W5lZGSc0v0BAAg214zJ0Mj0BFU3NOvhZTtMx7GFDped2tpaZWdny+VyHff2BQsWKD8/X7Nnz9bGjRuVnZ2tyZMn6/Dhw63b5OTkaMSIEd/6KS4ubt2mvLxc119/vZ5++ulT+LUAAAhOTkfLZGVJemXDARUUVZgNZAOWtxOLcViWpddee02XX35563Xjxo1Tbm6unnzySUmSx+NRRkaGbr31Vt19993tetyGhgZdcMEFmjlzpq677rqTbtvQ8OXXfFVVVcrIyFBlZaXi4+M7/ksBABAA8l8u0D82HlR2nwS99rOz5XBYpiP5VVVVlRISEvzy+e3TOTuNjY3asGGD8vLyvnwCh0N5eXlau3Ztux7D6/Xqhhtu0HnnnXfSoiNJ8+bNU0JCQusPu7wAAHZw95Rhio0M06YDlXp1wwHTcYKaT8tOWVmZ3G63UlJSvnZ9SkqKSkrat8DZ6tWrtWDBAi1cuFA5OTnKycnRxx9/3Ob299xzjyorK1t/ioqKOvU7AAAQCJLjonTb+YMlSb9dtkOVx5oMJwpeYaYDfNP48ePl8bT/3AKRkZGKjIz0YyIAAMyYflZ//d8H+7X7SK3+sPwTzZqaaTpSUPLpNztJSUlyOp0qLS392vWlpaVKTU315VMBAGB7EWEOzZnaMll5/tq92lVabThRcPJp2YmIiNDo0aO1YsWK1us8Ho9WrFihM88805dPBQBASDhnSC9dmJkit8erOYu2qhPHFYWsDpedmpoaFRQUqKCgQJK0Z88eFRQUaP/+/ZKk/Px8PfPMM5o/f762b9+um266SbW1tZoxY4ZPgwMAECp+dUmmIsIcWvPpUS3b0r45sPhSh+fsfPjhh5o0aVLr5fz8fEnS9OnT9fzzz+vaa6/VkSNHNGvWLJWUlCgnJ0fLli371qRlAADQPn17Ruun5wzUH98p1INLt2vi0GR1i3CajhU0OnWenUDicrnkcrnkdru1a9cuzrMDALCVY41unf/ouyqurNfPzx+s/AuGmI7kU/48z45tys4X/DlYAACYtHTzId380kZFhDm0Iv9cZSRGm47kM0FzUkEAAOA/F49M1ZkDe6qx2aMHl24zHSdoUHYAAAgSltWybpbTYemfW0v170+OmI4UFCg7AAAEkaGpcbrujH6SpDmLtqrJ3f4T8YYqyg4AAEHmvy4Yop4xEfr0SK3mr9lrOk7Ao+wAABBkErqF6xcXDZUkPb78Ex2urjecKLBRdgAACEJXj87QqD4Jqmlo1sPLdpqOE9BsU3ZcLpcyMzOVm5trOgoAAH7ncFiaOy1LYQ5LsZFh8nhsdSYZn+I8OwAABLFDlcfUO6Gb6Ridxnl2AADAcdmh6PgbZQcAANgaZQcAANgaZQcAANgaZQcAANgaZQcAANgaZQcAANiabcoOJxUEAADHw0kFAQCAcZxUEAAA4BRRdgAAgK1RdgAAgK1RdgAAgK1RdgAAgK2FmQ7ga18cXFZVVWU4CQAAaK8vPrf9cZC47crO0aNHJUkZGRmGkwAAgI46evSoEhISfPqYtis7iYmJkqT9+/f7dLByc3P1wQcf+Pw+bW3Tkeu/eV1bl6uqqpSRkaGioiKfnsPAH2Nzott9NTZf/XugjE1nXjMnuu1Ux8Zf43KirJ3Z3hdjY/rf04mydmb7jo5NIL7XtJWrs9vbYWy6+jOqrds6MzaVlZXq27dv6+e4L9mu7DgcLdOQEhISfPpCcjqdHX689tynrW06cv03rzvZ5fj4+IAfmxPd7quxOd72psemM6+ZE93W2bHx9bicKGtntvfF2Jj+93SirJ3ZvqNjE4jvNW3l6uz2dhibrv6Maus2X4zNF5/jvsQE5Xa6+eab/XKftrbpyPXfvO5kl33NH2Nzott9NTb+HpdTeY7OvGZOdBtj0/6xMf3v6VSewx9jE4jvNafyHKEyNl39GdXWbYE4NhLLRYQUxqZtjM3xMS5tY2zaxti0jbFpG8tFdEBkZKRmz56tyMhI01ECDmPTNsbm+BiXtjE2bWNs2sbYtM2fY2O7b3YAAAC+ynbf7AAAAHwVZQcAANgaZQcAANgaZQcAANgaZQcAANhayJadoqIiTZw4UZmZmRo1apReeeUV05ECyhVXXKEePXroqquuMh3FuCVLlmjo0KEaPHiwnn32WdNxAgqvk+Pj/aVtFRUVGjNmjHJycjRixAg988wzpiMFlLq6OvXr10933nmn6SgBpX///ho1apRycnI0adKkDt8/ZA89P3TokEpLS5WTk6OSkhKNHj1au3btUkxMjOloAeHdd99VdXW15s+fr1dffdV0HGOam5uVmZmplStXKiEhQaNHj9aaNWvUs2dP09ECAq+T4+P9pW1ut1sNDQ2Kjo5WbW2tRowYoQ8//JB/U5+77777VFhYqIyMDP3ud78zHSdg9O/fX1u2bFFsbOwp3T9kv9np3bu3cnJyJEmpqalKSkpSeXm52VABZOLEiYqLizMdw7j169crKytL6enpio2N1ZQpU/TWW2+ZjhUweJ0cH+8vbXM6nYqOjpYkNTQ0yOv1KkT/z/0tn3zyiXbs2KEpU6aYjmI7AVt2Vq1apalTpyotLU2WZWnhwoXf2sblcql///6KiorSuHHjtH79+lN6rg0bNsjtdisjI6OTqbtGV45NsOvsWBUXFys9Pb31cnp6ug4ePNgV0f2O11HbfDk2wfb+cjK+GJuKigplZ2erT58+uuuuu5SUlNRF6f3HF+Ny5513at68eV2UuOv4Ymwsy9K5556r3Nxcvfjiix3OELBlp7a2VtnZ2XK5XMe9fcGCBcrPz9fs2bO1ceNGZWdna/LkyTp8+HDrNl/sE/7mT3Fxces25eXluv766/X000/7/Xfyla4aGzvwxVjZFWPTNl+NTTC+v5yML8ame/fu2rRpk/bs2aOXXnpJpaWlXRXfbzo7Lq+//rqGDBmiIUOGdGXsLuGL18x7772nDRs2aNGiRXrooYe0efPmjoXwBgFJ3tdee+1r140dO9Z78803t152u93etLQ077x589r9uPX19d4JEyZ4X3jhBV9F7XL+Ghuv1+tduXKl97vf/a4vYgaEUxmr1atXey+//PLW22+77Tbviy++2CV5u1JnXkd2e51806mOjR3eX07GF+8/N910k/eVV17xZ8wudyrjcvfdd3v79Onj7devn7dnz57e+Ph479y5c7sydpfwxWvmzjvv9P75z3/u0PMG7Dc7J9LY2KgNGzYoLy+v9TqHw6G8vDytXbu2XY/h9Xp1ww036LzzztN1113nr6hdzhdjEyraM1Zjx47Vli1bdPDgQdXU1OjNN9/U5MmTTUXuMryO2taesbHr+8vJtGdsSktLVV1dLUmqrKzUqlWrNHToUCN5u0p7xmXevHkqKirS3r179bvf/U4zZ87UrFmzTEXuMu0Zm9ra2tbXTE1Njd555x1lZWV16HmCsuyUlZXJ7XYrJSXla9enpKSopKSkXY+xevVqLViwQAsXLlROTo5ycnL08ccf+yNul/LF2EhSXl6err76ar3xxhvq06ePLT/g2jNWYWFhevTRRzVp0iTl5OTojjvuCImjRtr7OgqF18k3tWds7Pr+cjLtGZt9+/ZpwoQJys7O1oQJE3Trrbdq5MiRJuJ2GV+9L9tRe8amtLRU48ePV3Z2ts444wxdf/31ys3N7dDzhPkscZAZP368PB6P6RgBa/ny5aYjBIxp06Zp2rRppmMEJF4nx8f7S9vGjh2rgoIC0zEC2g033GA6QkAZOHCgNm3a1KnHCMpvdpKSkuR0Or81qa20tFSpqamGUgUGxqb9GKu2MTZtY2zaxtgcH+PStq4am6AsOxERERo9erRWrFjRep3H49GKFSt05plnGkxmHmPTfoxV2xibtjE2bWNsjo9xaVtXjU3A7saqqalRYWFh6+U9e/aooKBAiYmJ6tu3r/Lz8zV9+nSNGTNGY8eO1eOPP67a2lrNmDHDYOquwdi0H2PVNsambYxN2xib42Nc2hYQY9Oxg8a6zsqVK72SvvUzffr01m2eeOIJb9++fb0RERHesWPHetetW2cucBdibNqPsWobY9M2xqZtjM3xMS5tC4SxCdm1sQAAQGgIyjk7AAAA7UXZAQAAtkbZAQAAtkbZAQAAtkbZAQAAtkbZAQAAtkbZAQAAtkbZAQAAtkbZAQAAtkbZAQAAtkbZAQAAtkbZAQAAtvb/AZ2vHqTjzpqJAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(hm.k, hm.tracer_profile_ukm[:, 1000])\n", "plt.xscale(\"log\")\n", "plt.yscale(\"log\")\n", "plt.xlim((1e-2, 1e5));\n", "# plt.ylim((1e-1,1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And we've set up our profile model.\n", "\n", "This profile model has 2 additional model parameters. You can always update these parameters using:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "hm.tracer_profile_params = {\"a\": 0.049, \"b\": 2.248}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Creating a new concentration-mass relation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The concentration-mass relation we use here follows the one from Maccio et al.(2007):" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "c_{\\rm HI}(M,z) = c_0 \\Big(\\frac{M}{10^{11}{\\rm M_\\odot h^{-1}}}\\Big)^{-0.109}\\frac{4}{(1+z)^\\gamma}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, because `halomod` already has a generic `CMRelation` class in place, all you really need is to specify a `cm` function for the equation above:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "from hmf.halos.mass_definitions import SOMean\n", "\n", "from halomod.concentration import CMRelation" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "class Maccio07(CMRelation):\n", " \"\"\"\n", " HI concentration-mass relation based on Maccio et al.(2007).\n", " Default value taken from 1611.06235.\n", " \"\"\"\n", "\n", " _defaults = {\"c_0\": 28.65, \"gamma\": 1.45}\n", " native_mdefs = (SOMean(),)\n", "\n", " def cm(self, m, z):\n", " return (\n", " self.params[\"c_0\"] * (m * 10 ** (-11)) ** (-0.109) * 4 / (1 + z) ** self.params[\"gamma\"]\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that for the concentration-mass relation that you put in, you need to specify the mass definition that this relation is defined in. In this case, it is defined in Spherical-Overdensity method, which is the default." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And set this to be the model for the tracer concentration-mass relation:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "hm.tracer_concentration_model = Maccio07" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And check the concentration-mass relation at z=0:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjYAAAGpCAYAAACebnw+AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAL6BJREFUeJzt3Xtw1HWC9/tPLiQhkARCSEJIuiPeuSWSpDvMODOiOIwzMorCgOmdZd05TpXPlLtVnJ2qcevs+uzz1DOeqtmzZc3aVVM7VTuuZQcQHVFxx5kVZViVdCCQACJEFLs7QBIC5Aq5df/OH006CfBTQi7d/cv7VZU/8gOTL/1D+22H8EkwDMMQAACABSRG+wAAAAAThbABAACWQdgAAADLIGwAAIBlEDYAAMAyCBsAAGAZhA0AALAMwgYAAFgGYQMAACyDsAEAAJZB2AAAAMuY8rBpb29XeXm5SktLtXTpUv32t7+d6iMAAACLSpjqEcxgMKi+vj6lp6erp6dHS5cu1YEDBzRv3rypPAYAALCgKX/FJikpSenp6ZKkvr4+GYYhBsYBAMBEGHPY7N27V2vXrlVBQYESEhK0c+fOa36O2+1WcXGx0tLS5HQ6VVtbO+rH29vbVVJSosLCQv385z9XTk7OTf8CAAAAhiSP9R/o6elRSUmJ/vqv/1qPPfbYNT++fft2bdmyRb/5zW/kdDr1wgsvaM2aNTpx4oRyc3MlSXPmzFFDQ4NaWlr02GOPaf369crLy7vu5+vr61NfX1/k/VAopAsXLmjevHlKSEgY6/EBAEAUGIahrq4uFRQUKDFxEr9gZIyDJOONN94Ydc3hcBg/+9nPIu8Hg0GjoKDAeP7556/7MZ5++mljx44dpp/jueeeMyTxxhtvvPHGG28WeAsEAuNJj6815ldsvkp/f7/q6ur07LPPRq4lJiZq9erV2rdvnySppaVF6enpysjIUEdHh/bu3aunn37a9GM+++yz2rJlS+T9jo4O2Ww2BQIBZWZmTuTxAQDAJOns7FRRUZEyMjIm9fNMaNi0tbUpGAxe82WlvLw8HT9+XJLk8/n005/+NPKHhp955hktW7bM9GOmpqYqNTX1muuZmZmEDQAAcWay/xjJhIbNjXA4HKqvr5/qTwsAAKaBCf3TOzk5OUpKSlJLS8uo6y0tLcrPz5/ITwUAAHCNCQ2blJQUlZWVaffu3ZFroVBIu3fv1sqVKyfyUwEAAFxjzF+K6u7u1smTJyPvnzp1SvX19crOzpbNZtOWLVu0efNmlZeXy+Fw6IUXXlBPT4+efPLJCT04AADA1cYcNgcOHNCqVasi7w99x9LmzZv10ksvaePGjTp37pz+8R//Uc3NzSotLdW7775r+vfUAAAATJQp34oar87OTmVlZamjo4PvigIAIE5M1fP3lG9F3Sy3263FixeroqIi2kcBAAAxildsAADApOMVGwAAgDEibAAAgGUQNgAAwDIIGwAAYBmEDQAAsAzCBgAAWAZhAwAALCNuwoa/oA8AAHwd/oI+AAAw6fgL+gAAAMaIsAEAAJZB2AAAAMsgbAAAgGUQNgAAwDIIGwAAYBmEDQAAsAzCBgAAWAZhAwAALCNuwoZJBQAA8HWYVAAAAJOOSQUAAIAxImwAAIBlEDYAAMAyCBsAAGAZhA0AALAMwgYAAFgGYQMAACyDsAEAAJZB2AAAAMsgbAAAgGXETdiwFQUAAL4OW1EAAGDSsRUFAAAwRoQNAACwDMIGAABYBmEDAAAsg7ABAACWQdgAAADLIGwAAIBlEDYAAMAyCBsAAGAZcRs2rx0IqKdvMNrHAAAAMSRuw+Z/vn1Mlb/crX9886iON3dG+zgAACAGJEf7ADfLlj1TTT2DenmfTy/v86nMPlcup03fX7ZAaTOSon08AAAQBXE7gnnxYruOnhuQx+vTfx1r0WAo/MuYkz5D61cUqspp06L5s6N8WgAAIE3dCGbchI3b7Zbb7VYwGFRjY+OoB6a1s1evHghoa21Ap9svR/6Zb9w6Ty6nXQ8uzlNKctx+1Q0AgLhH2Jj4qgcmGDK050SrPF6/PjjRqqFfWc7sVG2sKNSmCpuKstOjcGoAAKY3wsbEjT4wTRcvafv+gLbtD+hcV58kKSFBuu+O+XI57Vp1V66SEhOm6tgAAExrhI2JsT4wA8GQ3jvWIo/Xrw9PtkWuF2SlaZPDpo0VRcrLTJvMIwMAMO0RNibG88CcauvR1lq/dhwI6OKlAUlSUmKCHrw7T65Km755a44SeRUHAIAJR9iYmIgHpncgqHePNsvj9Wn/lxcj1+3z0lXlsGl9WaHmzU6dqCMDADDtETYmJvqBOdHcpWqvT78/eFpdV/4m45SkRD20LF8up10VxXOVkMCrOAAAjAdhY2KyHphL/YN6u+GMPF6/Djd1RK7fnjtbLqdN61YUKmvmjAn7fAAATCeEjYmpeGAON7Wr2uvXm/VndHkgKElKm5GoH5YUyOW0a3lhFq/iAAAwBoSNial6YCSps3dAOw+dlqfGrxMtXZHrSxdmqsph1yOlBZqVGrerFAAATBnCxsRUhs0QwzBU57soj9evd46cVf9gSJI0OzVZj94TfhXn7gVTcxYAAOIRYWMiGmEz0sWefr1W16TqWr9OtfVErq+wzZHLadcPljPCCQDA1QgbE9EOmyGhkKF9X5xXtdevP37SHBnhzJo5Q+vLwiOctzLCCQCAJMLGVKyEzUitXb3acaBJ1V7/qBHOlYvmyVVp03cX5zPCCQCY1ggbE7EYNkOCIUN7G8/J4/Xp/eOtCkVGOFP0o/IiPeFghBMAMD0RNiZiOWxGOt1+Wdtr/dq2P6DWESOc3xka4bxzvpKTeBUHADA9EDYm4iVshgwEQ9r9aXiE878/Gx7hXJCVpk0V4RHO/CxGOAEA1kbYmIi3sBnpyysjnK9eNcK5+u5cuZx23XsbI5wAAGsibK7idrvldrsVDAbV2NgYl2EzpHcgqD9+0ixPjV+1X16IXLdlp6vKadMGRjgBABZD2JiI51dsrqexpUvVXr9eP9ikrt7wCOeMpAQ9tHSBXE6bHLdkM98AAIh7hI0Jq4XNkEv9g9rVcFYer08NI0Y4b8udrSqHTY+vKFRWOiOcAID4RNiYsGrYjHT0dIc8Xp/erD+jS/3hEc7U5EStLSmQy2lTadEcXsUBAMQVwsbEdAibIZ29A3rz0Gl5vH4dbx4e4Vy8IFOuSpseKV2o2YxwAgDiAGFjYjqFzRDDMHTQHx7h3HV4eIRzVkqSHr1noVxOuxYXTI/HAgAQnwgbE9MxbEa62NOv1w+G5xu+GDHCec+VEc6HGeEEAMQgwsbEdA+bIYYRHuH0eP3649HRI5yPrwiPcN6WywgnACA2EDYmCJtrmY1wVi7Klstp15oljHACAKKLsDFB2JgLhgzt/eycPDV+vX+8ZdQI54byIlUxwgkAiBLCxgRhc2POtF/Wtv0Bbav1jxrh/Pbt8+Vy2nT/XbmMcAIApgxhY4KwGZvwCGerPF7fqBHO/Mw0bXIUaVOFjRFOAMCkI2xMEDY3z3e+R9W1fu040KQLPf2SwiOcD9yVK1elXd9ihBMAMEkIGxOEzfj1DQb17tFmebx+1Z4aHuEsyp6pKoddG8oLlcMIJwBgAhE2JgibifVZS5c81xnhXLMkX39RaZeTEU4AwAQgbEwQNpPjcn9Qbx8+I4/Xr4ZAe+T6rfNnyeW0M8IJABgXwsYEYTP5wiOcfr1Zf3rUCOfDywvkqrTpHkY4AQBjRNiYIGymTlfvgN6sP6NXanyjRjjvXpApl9OmR+9hhBMAcGMIGxOEzdQzDEOHAu3y1Pi16/AZ9Y0Y4XzknoVyOW1aUpAV5VMCAGIZYWOCsImu9kv9ev3gaXm8Pn1xbniEs7RojlxOmx5eXqCZKYxwAgBGI2xMEDaxwTAM1XxxQR6vT3/8pFkDwfBvo8y0ZD1eViiX06bbcjOifEoAQKwgbEwQNrHnXFefdtQFVO31q+ni8Ain85ZsuSrtWrMkT6nJvIoDANMZYWOCsIldoaERTq9fuz8dHuGcN2t4hNM2jxFOAJiOCJuruN1uud1uBYNBNTY2EjYx7mzHZW2rDWjbfr9aOvsi1799R3iE8wFGOAFgWiFsTPCKTXwZDIa0+3irPF6/9jaei1zPy0zVpgqbNjmKtCBrZhRPCACYCoSNCcImfvnPX7oywhnQ+SsjnIkJ0gN358nltOnbt89nhBMALIqwMUHYxL++waD++EmLPDU+eUeMcBbOnakqp00byoo0P4MRTgCwEsLGBGFjLSdbr4xw1jWpc8QI53eX5MvltGnlonnMNwCABRA2Jggba7rcH9Q7R87K4/XpkL89cn3R/Fmqcti0vqxQc9JTondAAMC4EDYmCBvr++RMh6q9fu08dFo9I0Y4f7B8gVxOu1bYGOEEgHhD2JggbKaP7r5BvVl/Wq/U+PXp2c7I9bvyM+SqtOvR0gJlpM2I4gkBADeKsDFB2Ew/hmGoPtAuj9evtxuGRzjTU5L0SGl4hHPpQkY4ASCWETYmCJvprePSgF4/2CSP16fPR4xwllwZ4VzLCCcAxCTCxgRhAyn8Ko731AV5vH69e/RsZIQzIy1Zj68Ij3DenscIJwDECsLGBGGDq7V192nHgSZV1/oUuDA8wum4JVsup03fW5rPCCcARBlhY4KwgZlQyNB/n2yTp8an3cdbFbyywpk9K0UbygtV5bDJPm9WlE8JANMTYWOCsMGNONtxWdv3B7StNqDmzt7I9W/dniOX064H7s7VDEY4AWDKEDYmCBuMxWAwpPeHRjg/O6eh3+15manaWGHTpooiFcxhhBMAJhthY4Kwwc0KXLikrbV+vXogoLbu4RHO+++6MsJ5x3wlMcIJAJOCsDFB2GC8+gdD+tOxZnlq/Nr3xfnI9YVzroxwlhcqNyMtiicEAOshbEwQNphIJ1u7tbXWr9fqmtRxeUCSlJyYoDVDI5y3MsIJABOBsDFB2GAy9A4E9c7h8AjnwZEjnDmzVOW06fEVhZo7ixFOALhZhI0JwgaT7diZTlXX+vTGweERzpTkRD28bIFclTatsM3lVRwAGCPCxgRhg6nS3Teot+rP6JUan45dPcLptOnRexYywgkAN4iwMUHYYKoZhqGGpg55anx6+/AZ9Q6MHOEskMtpZ4QTAL4GYWOCsEE0dVwa0O8PNcnj9etka3fkeklhllxOux4uWaD0lOQonhAAYhNhY4KwQSwwDEO1V0Y4/3CdEc4qp013MMIJABGEjQnCBrGmrbtPr9U1qdrrl//Cpch1R3G2XJWMcAKARNiYImwQq0IhQx+ebJPH69N7n141wllWqCccNhXnMMIJYHoibEwQNogHLZ292r4/oK21fp3tGD3CWeWwafXiPEY4AUwrhI0JwgbxZDAY0gcnzsnj9enPjcMjnPMzUrWpokibHDYtZIQTwDRA2JggbBCvAhcuadt+v7bvv3qEM1cup50RTgCWRtiYIGwQ7/oHQ/qvYy3yeH36+PPRI5xPOIr0o4oiRjgBWA5hY4KwgZV8fq5bW71+7bhqhPO7S/Lkctq1ctE8JfIqDgALIGxMEDawot6BoP7zyFl5vH7V+S5Grt+SM0tVDpvWlzHCCSC+ETZXcbvdcrvdCgaDamxsJGxgWZ+e7VS11683Dp1Wd9+gpPAI5w+WLZDLaVOZnRFOAPGHsDHBKzaYLnr6BvVWQ3iE85MzwyOcd+ZlyFUZHuHMZIQTQJwgbEwQNphuDMPQ4aYOebw+vdUwPMI5c8bwCOeyQkY4AcQ2wsYEYYPprOPygN44GB7h/GzECOfywiy5nDatLSlghBNATCJsTBA2QPhVnP1fXpTH69MfjjSrPxh+FScjNVmPrVioKqddd+YzwgkgdhA2JggbYLTzQyOctX75zg+PcFYUz5XLadf3luYrbQYjnACii7AxQdgA1xcKGfro8zZVe/3607GWyAjn3PQZ2lBepCccNt3CCCeAKCFsTBA2wNdr6ezVq1dGOM+MGOG897YcuZyMcAKYeoSNCcIGuHHBkKE9J1rl8fr1wYnWUSOcG8uLtMlRpMK56dE9JIBpgbAxQdgANydw4ZK27w9o2/6A2rr7JEkJCdKqO3Plctp03525jHACmDSEjQnCBhifgeDwCOdHJ0ePcG6qKNLGiiLlZjLCCWBiETYmCBtg4nxxrltba8MjnO2Xhkc4H1wcHuH8xq2McAKYGISNCcIGmHi9A0H94ehZeWr8OjBihLN4XrqqnDatLytSNiOcAMaBsDFB2ACT63hzeITz9wdHjHAmJer7y/LlqrSrnBFOADeBsDFB2ABTo6dvUG83nJHH69eR0x2R63fkzZbLade6FYxwArhxhI0JwgaYeoeb2uWp8euthjO6PBCUFB7h/GFJgVyVNi0vnBPdAwKIeYSNCcIGiJ6OywPaeei0PF6fGluGRziXLQyPcP6wlBFOANdH2JggbIDoMwxDB3wXVe31653DZ0eNcK5bsVBVTpvuyuffTwDDCBsThA0QWy709Ov1uiZ5vD59OWKEs9w+V65Kmx5auoARTgCEjRnCBohNoZChfV+cl8fr058+adHglRHOOekztH5FoaqcNi2aPzvKpwQQLYSNCcIGiH2tnb169UBAW2sDOt1+OXL9G7fOk8tp14OL85SSzAgnMJ0QNiYIGyB+BEOG/tzYKk+NX++PGOHMmZ2qjRWF2lRhU1E2I5zAdEDYmCBsgPjUdHF4hPNc1/AI5313zJfLadequxjhBKyMsDFB2ADxbSAY0nvHWuTx+vXhybbI9YKsNG1y2LSxokh5jHAClkPYmCBsAOs41dYTHuE8ENDFKyOcSYkJevDuPLkqbfrmrTmMcAIWQdiYIGwA6+kdCOrdo83yeH3a/+XwCKd9XrqqHDatLyvUvNmpUTwhgPEibEwQNoC1nWjuUrXXp98fPK2uESOcDy3Ll8tpV0UxI5xAPCJsTBA2wPRwqX94hPNw0/AI5+25s+Vy2rRuRaGyZjLCCcQLwsYEYQNMP4eb2lXt9evN+uERzrQZieERTqddywuzeBUHiHGEjQnCBpi+OnsH9Oah03qlxq8TLV2R60sXZsrltOuHJQWalcoIJxCLCBsThA0AwzB00H9Rnhq/dh05q/7B8Ajn7NRkrbsnPMJ59wL++wDEEsLGBGEDYKSLPf16/WCTPF6/TrX1RK6vsM2Ry2nXD5YzwgnEAsLGBGED4HoMw9C+z8/L4/Xrj580R0Y4s2bO0Pqy8AjnrYxwAlFD2JggbAB8ndauXu040KRqr3/UCOfKRfPkqrTpu4vzGeEEphhhY4KwAXCjgiFDexvPyeP16f3jrQpFRjhT9KPyIj3hYIQTmCqEjQnCBsDNON1+Wdtr/dq2P6DWESOc3xka4bxzvpKTeBUHmCyEjQnCBsB4DARD2v1peITzvz8bHuFckJWmTRXhEc78LEY4gYlG2JggbABMlC+vjHC+etUI5+q7c+Vy2nXvbYxwAhOFsDFB2ACYaH2DQyOcftWeuhC5bstOV5XTpg2McALjRtiYIGwATKbPWrrk8fr1+sEmdfUOj3B+b2m+XE6bHLdkM98A3ATCxgRhA2AqXO4P6u3D4RHOhkB75PptV0Y4H7unUFnpjHACN4qwMUHYAJhqR093yOP1683607rUPzzCuXZ5gaqcNpUWzeFVHOBrEDYmCBsA0dLVO6Cd9WfkqfHpePPwCOfiBZlyVdr0SOlCzWaEE7guwsYEYQMg2sIjnO3yeH3adXj0COej9xSoymHX4gL++wSMRNiYIGwAxJKhEc5qr19fjBjhvOfKCOfDjHACkggbU4QNgFhkGIb2fXFlhPPo6BHOx1eERzhvy2WEE9MXYWOCsAEQ64ZGOLfW+tV0cXiEs3JRtlxOu9YsYYQT0w9hY4KwARAvgiFDez87J0+NX+8fbxk1wrmhvEhVjHBiGiFsTBA2AOLRmfbL2rY/oG21/lEjnN++fb5cTpvuvyuXEU5YGmFjgrABEM/CI5yt8nh9o0Y48zPTtMlRpE0VNkY4YUmWDZtAIKAf//jHam1tVXJysv7hH/5BGzZsuOF/nrABYBW+8z3aWhvQqwcCutDTLyk8wvnAXblyVdr1LUY4YSGWDZuzZ8+qpaVFpaWlam5uVllZmRobGzVr1qwb+ucJGwBWYzbCWZQ9U1UOuzaUFyqHEU7EOcuGzdVKSkq0a9cuFRUV3dDPJ2wAWNlnLV2qrvXr9bomdV4Z4ZyRlKDvLV0gl9MmJyOciFNT9fw95j+ptnfvXq1du1YFBQVKSEjQzp07r/k5brdbxcXFSktLk9PpVG1t7XU/Vl1dnYLB4A1HDQBY3e15GXpu7RJ5/361frV+uUqL5mggaOjthjPa9G81Wv0vf9a/f3hKHZcGon1UICaNOWx6enpUUlIit9t93R/fvn27tmzZoueee04HDx5USUmJ1qxZo9bW1lE/78KFC/rLv/xL/du//dvNnRwALGxmSpI2lBdp58++qV3P3Ksqp03pKUn6/FyP/teuY3L88j393Y4GHfJfVJx9Dwgwqcb1paiEhAS98cYbevTRRyPXnE6nKioq9OKLL0qSQqGQioqK9Mwzz+gXv/iFJKmvr08PPvignnrqKf34xz/+ys/R19envr6+yPudnZ0qKiriS1EApp2u3gG9WX9Gr1xnhLPKadOj9zDCidgVs1+K+ir9/f2qq6vT6tWrhz9BYqJWr16tffv2SQr/teN/9Vd/pfvvv/9ro0aSnn/+eWVlZUXe+LIVgOkqI22G/qLSrj/87bf0+//xDT2+olCpyYk6drZT/8/Oo3L+n/f0928c0SdnOqJ9VCBqJjRs2traFAwGlZeXN+p6Xl6empubJUkfffSRtm/frp07d6q0tFSlpaU6cuSI6cd89tln1dHREXkLBAITeWQAiDsJCQlaYZur/+9HJfL+/QP6h4cXa9H8WerpD6ra69cPfv2hHnV/pB0HArrcH4z2cYEpNeWvWd57770KhUI3/PNTU1OVmsq3OQLA9cxJT9FP7r1Ff/3NYtV8cUEer09//KRZ9YF21Qfa9b93HdPjZYVyOW26LTcj2scFJt2Ehk1OTo6SkpLU0tIy6npLS4vy8/Mn8lMBAEZISEjQylvnaeWt83Suq0876gKq9oZHOH/30Zf63UdfynlLtlyVdq1ZkqfU5KRoHxmYFBP6paiUlBSVlZVp9+7dkWuhUEi7d+/WypUrJ/JTAQBMzM9I1f+47zbt/fkqvfRkhR5cnKfEBMl76oL+ZushfeP59/X//uG4/OcvRfuowIQb8ys23d3dOnnyZOT9U6dOqb6+XtnZ2bLZbNqyZYs2b96s8vJyORwOvfDCC+rp6dGTTz45oQcHAHy1xMQE3Xdnru67M1dnOy5rW21A2/b71dLZp9/8+XP95s+f69t3hEc4H2CEExYx5m/33rNnj1atWnXN9c2bN+ull16SJL344ov61a9+pebmZpWWlurXv/61nE7nhByYv3kYAG7eYDCk3cdb5fH6tbfxXOR6XmaqNlXYtMlRpAVZM6N4QljVtJlUGCvCBgAmhv/8JW3d79er+wM6f2WEMzFBeuDuPLmcNn379vmMcGLCEDZXcbvdcrvdCgaDamxsJGwAYIL0DQb1p09a5PH6VPPF6BHOJxw2bSgr0vwMvjsV40PYmOAVGwCYPCdbu1TtDei1usCoEc41S/LlctpVuYgRTtwcwsYEYQMAk693IKhdh8/K4/XpkL89cn3R/Fmqcti0vqxQc9JTondAxB3CxgRhAwBT65MzHar2+rXz0Gn1XPmbjFOTE/WD5Qvkctq1wjaHV3HwtQgbE4QNAERHd9+g3qw/rVdq/Pr0bGfk+l35GXJV2vVoaYEy0mZE8YSIZYSNCcIGAKLLMAzVB9rl8fr1dsMZ9Q2GZ3LSU5L0SOlCuZw2LV2YFeVTItYQNiYIGwCIHR2XBvT6wSZ5vD59fq4ncr2kaI5cTpvWLi/QzBTmG0DYmCJsACD2GIYh76kL8nj9evfoWQ0Ew08tGWnJenxFeITz9jxGOKczwsYEYQMAsa2tu087DjSputanwIXLkeuOW7Llctr0vaX5jHBOQ4SNCcIGAOJDKGTov0+2qdrr03uftioYCj/dZM9K0YbyQlU5bLLPmxXlU2KqEDZX4W8eBoD41dzRq+37A9pa61dzZ2/k+rduz5HLadfquxnhtDrCxgSv2ABA/BoMhvTBiXPyeH36c+M5DT0D5WWmamOFTZsqilQwhxFOKyJsTBA2AGANgQuXtLXWr1cPBNTWPTzCef9deXJVhkc4kxjhtAzCxgRhAwDW0j8Y0p+ONctT49e+L85HrhfODY9w/qicEU4rIGxMEDYAYF0nW7u1tdav1+qa1HF5QJKUnDg0wmnTylvnMd8QpwgbE4QNAFhf70BQ71wZ4Tw4coQzZ5aqnIxwxiPCxgRhAwDTy7Eznaqu9emNg8MjnCnJiXp42QK5Km1aYZvLqzhxgLAxQdgAwPTU3Teot+rPyOP16ZMzV41wOm169J6FjHDGMMLGBGEDANObYRhqaOqQp8antw+fUe/AyBHOArmcdkY4YxBhY4KwAQAM6bg0oN8fapLH69fJ1u7I9ZLCLLmcdj1cskDpKclRPCGGEDYmCBsAwNUMw1BtZISzWf3B8Ks4QyOcVU6b7mCEM6oIm6swqQAAuBHnu/v0Wl2Tqmv98p2/FLnuKM6Wq5IRzmghbEzwig0A4EaEQoY++rxNnhq//uvTltEjnGWFesJhU3EOI5xThbAxQdgAAMaqpXN4hPNsx9UjnDY9cHeeZjDCOakIGxOEDQDgZg0GQ9pzZYRzz4gRztyMVG2qKNJGh00LGeGcFISNCcIGADARAhcuadt+v7bvb1Jbd5+koRHOXLmcdn37DkY4JxJhY4KwAQBMpP7BkP7rWIs8Xp8+/nx4hHPhnJmqctq0obxQuRlpUTyhNRA2JggbAMBk+fxct7Z6/dpx1Qjnd5fkyeW0a+WieUrkVZybQtiYIGwAAJOtdyCo/zxyVh6vX3W+i5Hrt+TMUpUjPMI5dxYjnGNB2JggbAAAU+nTs52q9vr1xqHT6u4blBQe4fzBsgVyOW0qszPCeSMIGxOEDQAgGnr6BvVWwxm9UjN6hPPOvAy5KsMjnJmMcJoibEwQNgCAaDIMQ4ebOuTx+vRWw/AI58wZwyOcywoZ4bwaYWOCsAEAxIqOywPaeei0Xqnx6bMRI5zLC7Pkctq0tqSAEc4rCJursBUFAIhVhmHogO+iPDU+/eeRESOcqcl6bMVCVTntujN/eo9wEjYmeMUGABDLLvT067W6gDze0SOcFcVz5XLa9b2l+UqbMf1GOAkbE4QNACAehEKGPv78vDxen/50bHiEc276DG0oL9ITDptumUYjnISNCcIGABBvWjp79eqVEc4zI0Y4770tPMK5erH1RzgJGxOEDQAgXgVDhvacaJXH69cHJ1ojI5zzr4xwbrLwCCdhY4KwAQBYQdPFS9pWG9C2/YFRI5yr7sxVldOm++7MtdQIJ2FjgrABAFjJQHB4hPOjk6NHODdVFGljRZFyM+N/hJOwMUHYAACs6otz3dpaGx7hbL80PML54OLwCOc3bo3fEU7CxgRhAwCwut6BoP5w9Kw8NX4dGDHCWTwvXVVOm9aXFSk7zkY4CRsThA0AYDo50dylaq9Pvz94Wl1DI5xJifr+sny5Ku0qj5MRTsLGBGEDAJiOLvUP6u2GM3qlxq8jpzsi1+/Imy2X0651K2J7hJOwMUHYAACmu8NN7ar2+vVm/RldHghKCo9w/rCkQK5Km5YXzonuAa+DsDFB2AAAENbZOzzC2dgyPMK5bGF4hPOHpbEzwknYmCBsAAAYzTAM1fkuyuP1653DZ0eNcK5bsVBVTpvuyo/ucyZhY4KwAQDA3IWefr1e1ySP16cvR4xwltvnylVp00NLF0RlhJOwuYrb7Zbb7VYwGFRjYyNhAwDAVwiFDO374soI5yctGrwywjknfYY2lBXqCYdNi+bPnrLzEDYmeMUGAICxae3s1asHAtpaG9Dp9suR69+8bZ6qHHY9uDhPKcmTO8JJ2JggbAAAuDnBkKE/N7bKU+PX+yNGOHNmp2pjRaE2VdhUlJ0+KZ+bsDFB2AAAMH5NFy9p+/7wCOe5rvAIZ0KCdN8d8/UXlfYJH+EkbEwQNgAATJyBYEjvHWuRx+vXhyfbItcLstK0yWHTxooi5U3ACCdhY4KwAQBgcnzZ1qOttX69eiCgi1dGOJMSE/Tg3XlyVdr0zVtzbnqEk7AxQdgAADC5egeCevdoszxen/Z/OTzCaZ+XriqHTevLCjVvduqYPiZhY4KwAQBg6jS2dKna69frdU2jRjgfWpYvl9OuiuIbG+EkbEwQNgAATL1L/YPa1XBWHq9PDU3DI5y3586Wy2nTuhWFypppPsJJ2JggbAAAiK4jTR2qrvVp56HhEc60GYnhEU6nXcsLs655FYewMUHYAAAQGzp7B/TmodN6pcavEy1dketLF2bK5bTrhyUFmpUaHuEkbEwQNgAAxBbDMHTQf1GeGr92HTmr/sHwCOfs1GStuyc8wrlwlgib6yFsAACIXRd7+vX6wSZ5vH6dauuJXF+em6K3/+/vTvrzd/KkfWQAADDtzJ2Vov/rW4v0k3tv0b7Pz8vj9euPnzSrPtA+JZ+fsAEAABMuISFB37gtR9+4LUetXb16+c+f6ucvTP7nndwpTwAAMO3lZqTpp9++dUo+F2EDAAAsI27Cxu12a/HixaqoqIj2UQAAQIziu6IAAMCkm6rn77h5xQYAAODrEDYAAMAyCBsAAGAZhA0AALAMwgYAAFgGYQMAACyDsAEAAJZB2AAAAMsgbAAAgGUQNgAAwDIIGwAAYBmEDQAAsAzCBgAAWAZhAwAALIOwAQAAlkHYAAAAyyBsAACAZRA2AADAMggbAABgGYQNAACwDMIGAABYBmEDAAAsI27Cxu12a/HixaqoqIj2UQAAQIxKMAzDiPYhxqKzs1NZWVnq6OhQZmZmtI8DAABuwFQ9f8fNKzYAAABfh7ABAACWQdgAAADLIGwAAIBlEDYAAMAyCBsAAGAZhA0AALAMwgYAAFgGYQMAACyDsAEAAJZB2AAAAMsgbAAAgGUQNgAAwDIIGwAAYBmEDQAAsAzCBgAAWAZhAwAALIOwAQAAlkHYAAAAyyBsAACAZRA2AADAMggbAABgGYQNAACwDMIGAABYBmEDAAAsg7ABAACWQdgAAADLIGwAAIBlEDYAAMAyCBsAAGAZcRM2brdbixcvVkVFRbSPAgAAYlSCYRhGtA8xFp2dncrKylJHR4cyMzOjfRwAAHADpur5O25esQEAAPg6hA0AALAMwgYAAFgGYQMAACyDsAEAAJZB2AAAAMsgbAAAgGUQNgAAwDIIGwAAYBmEDQAAsAzCBgAAWAZhAwAALIOwAQAAlkHYAAAAyyBsAACAZRA2AADAMggbAABgGYQNAACwDMIGAABYBmEDAAAsg7ABAACWQdgAAADLIGwAAIBlEDYAAMAyCBsAAGAZhA0AALAMwgYAAFgGYQMAACyDsAEAAJZB2AAAAMsgbAAAgGUQNgAAwDIIGwAAYBmEDQAAsAzCBgAAWAZhAwAALIOwAQAAlkHYAAAAyyBsAACAZRA2AADAMggbAABgGYQNAACwDMIGAABYBmEDAAAsg7ABAACWQdgAAADLIGwAAIBlEDYAAMAyohI269at09y5c7V+/fpofHoAAGBRUQmbv/3bv9XLL78cjU8NAAAsLCphc9999ykjIyManxoAAFjYmMNm7969Wrt2rQoKCpSQkKCdO3de83PcbreKi4uVlpYmp9Op2traiTgrAADAVxpz2PT09KikpERut/u6P759+3Zt2bJFzz33nA4ePKiSkhKtWbNGra2t4z4sAADAV0ke6z/w0EMP6aGHHjL98X/5l3/RU089pSeffFKS9Jvf/EbvvPOO/v3f/12/+MUvxnzAvr4+9fX1Rd7v6OiQJHV2do75YwEAgOgYet42DGNSP8+Yw+ar9Pf3q66uTs8++2zkWmJiolavXq19+/bd1Md8/vnn9U//9E/XXC8qKrrpcwIAgOg4f/68srKyJu3jT2jYtLW1KRgMKi8vb9T1vLw8HT9+PPL+6tWr1dDQoJ6eHhUWFmrHjh1auXLldT/ms88+qy1btkTeb29vl91ul9/vn9QHBl+vs7NTRUVFCgQCyszMjPZxpjXuRezgXsQW7kfs6OjokM1mU3Z29qR+ngkNmxv13nvv3fDPTU1NVWpq6jXXs7Ky+E0aIzIzM7kXMYJ7ETu4F7GF+xE7EhMn9xuyJ/Sj5+TkKCkpSS0tLaOut7S0KD8/fyI/FQAAwDUmNGxSUlJUVlam3bt3R66FQiHt3r3b9EtNAAAAE2XMX4rq7u7WyZMnI++fOnVK9fX1ys7Ols1m05YtW7R582aVl5fL4XDohRdeUE9PT+S7pMYrNTVVzz333HW/PIWpxb2IHdyL2MG9iC3cj9gxVfciwRjj913t2bNHq1atuub65s2b9dJLL0mSXnzxRf3qV79Sc3OzSktL9etf/1pOp3NCDgwAAGBmzGEDAAAQq6KyFQUAADAZCBsAAGAZhA0AALAMwgYAAFiGpcKmuLhYy5cvV2lp6XW/cwtT58SJEyotLY28zZw5Uzt37oz2saatf/7nf9aSJUu0dOlSvfLKK9E+zrSzbt06zZ07V+vXr7+h65g813vM29vbVV5ertLSUi1dulS//e1vo3jC6cPs9/94n8st9V1RxcXFOnr0qGbPnh3to2CE7u5uFRcXy+fzadasWdE+zrRz5MgRbd68WR9//LEMw9CqVav07rvvas6cOdE+2rSxZ88edXV16T/+4z/02muvfe11TJ7rPebBYFB9fX1KT09XT0+Pli5dqgMHDmjevHlRPq21mf3+H+9zuaVesUFseuutt/TAAw8QNVHy6aefauXKlUpLS9PMmTNVUlKid999N9rHmlbuu+8+ZWRk3PB1TJ7rPeZJSUlKT0+XJPX19ckwDFno//lj1mT9/o+ZsNm7d6/Wrl2rgoICJSQkXPfLFm63W8XFxUpLS5PT6VRtbe2oH09ISNB3vvMdVVRUyOPxTNHJrWki7seQV199VRs3bpzkE1vXeO/F0qVLtWfPHrW3t+vixYvas2ePTp8+PYW/gvg2kf8uYHwm8160t7erpKREhYWF+vnPf66cnJwJPr21TOa9GO9zecyETU9Pj0pKSuR2u6/749u3b9eWLVv03HPP6eDBgyopKdGaNWvU2toa+Tkffvih6urq9NZbb+mXv/ylDh8+PFXHt5yJuB+S1NnZqY8//ljf//73p+LYljTee7F48WL9zd/8je6//3499thjqqysVFJS0lT+EuLaRP27gPGbzHsxZ84cNTQ06NSpU6qurr5mzBmjTea9GPdzuRGDJBlvvPHGqGsOh8P42c9+Fnk/GAwaBQUFxvPPP3/dj/F3f/d3xu9+97tJPOX0MZ778fLLLxsul2sqjjktTMS/Gz/5yU+MXbt2TeYxLWs8j/8HH3xgPP7449d8TLPr+GqTcS+GPP3008aOHTsm9LxWNpn34maey2PmFZuv0t/fr7q6Oq1evTpyLTExUatXr9a+ffskheuxq6tLUvgPq77//vtasmRJVM5rdTdyP4bwZajJdaP3Yuj/kk6cOKHa2lqtWbNmys9qRWP5dwGTazz3oqWlJfL80dHRob179+rOO++c1PNa2XjuxUQ8l4953Tsa2traFAwGlZeXN+p6Xl6ejh8/Lin8G3PdunWSwn/C/amnnlJFRcWUn3U6uJH7IYX/A1FbW6vXX399qo84bdzovXjkkUfU0dGhWbNm6Xe/+52Sk+PiX/2Yd6OP/+rVq9XQ0KCenh4VFhZqx44dWrlypel1jN147kVSUpJ++tOfRv7Q8DPPPKNly5ZN9S/BMsZzL/Ly8sb9XG6Z/7otWrRIDQ0N0T4GRsjKyuLr1DGCVw+i67333hvTdUwes8e8vr5+ag8C03sx3ufyuPhSVE5OjpKSkq55kmxpaVF+fn6UTjV9cT9iB/ciunj8Ywf3InZE+17ERdikpKSorKxMu3fvjlwLhULavXs3L9tGAfcjdnAvoovHP3ZwL2JHtO9FzHwpqru7WydPnoy8f+rUKdXX1ys7O1s2m01btmzR5s2bVV5eLofDoRdeeEE9PT168skno3hq6+J+xA7uRXTx+McO7kXsiOl7MabvoZpEH3zwgSHpmrfNmzdHfs6//uu/GjabzUhJSTEcDodRU1MTvQNbHPcjdnAvoovHP3ZwL2JHLN8LS21FAQCA6S0u/owNAADAjSBsAACAZRA2AADAMggbAABgGYQNAACwDMIGAABYBmEDAAAsg7ABAACWQdgAAADLIGwAAIBlEDYAAMAyCBsAAGAZ/z+RQeiVKhh7KQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(hm.m, hm.tracer_concentration.cm(hm.m, 0))\n", "plt.xscale(\"log\")\n", "plt.xlim(1e5, 1e15)\n", "plt.ylim((1e1, 1e3))\n", "plt.yscale(\"log\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that this model has two additional parameters, which can be updated using:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "hm.tracer_concentration_params = {\"c_0\": 28.65, \"gamma\": 1.45}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "See the documentation for [`concentration.py`](https://halomod.readthedocs.io/en/docs/_autosummary/halomod.concentration.html) for more." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Creating a new HOD" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The HI HOD we use here is:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{split}\n", "\\langle M_{\\rm HI}^{\\rm cen}(M_h) \\rangle = M_h& \\Bigg[a_1^{\\rm cen}\\bigg(\\frac{M_h}{10^{10} M_\\odot}\\bigg)^{\\beta_{\\rm cen}} {\\rm exp}\\Big[{-\\bigg(\\frac{M_h}{M^{\\rm cen}_{\\rm break}}\\bigg)^{\\alpha_{\\rm cen}}}\\Big] \\\\&+a_2^{\\rm cen}\\Bigg] {\\rm exp}\\Big[{-\\bigg(\\frac{M_{\\rm min}^{\\rm cen}}{M_h}\\bigg)^{0.5}}\\Big]\n", "\\end{split}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{split}\n", "\\langle M_{\\rm HI}^{\\rm sat}(M_h) \\rangle = \n", "M_0^{\\rm sat}\\bigg( \\frac{M_h}{M^{\\rm sat}_{\\rm min}}\\bigg)^{\\beta_{\\rm sat}}\n", "{\\rm exp}\\Big[{-\\bigg(\\frac{M^{\\rm sat}_{\\rm min}}{M_h}\\bigg)^{\\alpha_{\\rm sat}}}\\Big]\n", "\\end{split}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the HOD, it's a bit more complicated. First, one needs to decide what type of tracer it is. The most generic class to use is `HOD`, however, you may prefer `HODBulk` where the tracer is considered to be continuously distributed ; or `HODPoisson`, which assumes Poisson distributed discrete satellite components, which is commonly used for galaxies.\n", "\n", "Second, if your model has a minimum halo mass to host any tracer as a sharp cut-off (or not), you should specify in your model definition\n", "\n", "```python\n", "sharp_cut = True # or False\n", "```\n", "\n", "If your model has a seperation of central and satellite components, you need to specify if the satellite occupation is inherently dependent on the existence of central galaxies:\n", "\n", "```python\n", "central_condition_inherent = False # or True\n", "```\n", "\n", "If False, the actual satellite component will be your satellite occupation times occupation for central galaxies.\n", "\n", "See the documentation for [`hod.py`](https://halomod.readthedocs.io/en/docs/_autosummary/halomod.hod.html), or the [reference paper](https://arxiv.org/abs/2009.14066) for more." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, you need to specify how to convert units between your HOD and the resulting power spectrum. For example, for HI the HOD is written in mass units, whereas the power spectrum is in temperature units. This is done by specifying a `unit_conversion` method." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Additionally, sometimes your HOD contains methods that need to be calculated, such as virial velocity of the halos, which you can just put into the class." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "import astropy.constants as astroconst\n", "import scipy.constants as const\n", "\n", "from halomod.hod import HODPoisson" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And in our case for utility, we also specify the number of galaxies in this model, which is defined by `_tracer_per_central` and `_tracer_per_satellite`:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "class Spinelli19(HODPoisson):\n", " \"\"\"\n", " Six-parameter model of Spinelli et al. (2019)\n", " Default is taken to be z=1(need to set it up manually via hm.update)\n", " \"\"\"\n", "\n", " _defaults = {\n", " \"a1\": 0.0016, # gives HI mass amplitude of the power law\n", " \"a2\": 0.00011, # gives HI mass amplitude of the power law\n", " \"alpha\": 0.56, # slope of exponential break\n", " \"beta\": 0.43, # slope of mass\n", " \"M_min\": 9, # Truncation Mass\n", " \"M_break\": 11.86, # Characteristic Mass\n", " \"M_1\": -2.99, # mass of exponential cutoff\n", " \"sigma_A\": 0, # The (constant) standard deviation of the tracer\n", " \"M_max\": 18, # Truncation mass\n", " \"M_0\": 8.31, # Amplitude of satellite HOD\n", " \"M_break_sat\": 11.4, # characteristic mass for satellite HOD\n", " \"alpha_sat\": 0.84, # slope of exponential cut-off for satellite\n", " \"beta_sat\": 1.10, # slope of mass for satellite\n", " \"M_1_counts\": 12.851,\n", " \"alpha_counts\": 1.049,\n", " \"M_min_counts\": 11, # Truncation Mass\n", " \"M_max_counts\": 15, # Truncation Mass\n", " \"a\": 0.049,\n", " \"b\": 2.248,\n", " \"eta\": 1.0,\n", " }\n", " sharp_cut = False\n", " central_condition_inherent = False\n", "\n", " def _central_occupation(self, m):\n", " alpha = self.params[\"alpha\"]\n", " beta = self.params[\"beta\"]\n", " m_1 = 10 ** self.params[\"M_1\"]\n", " a1 = self.params[\"a1\"]\n", " a2 = self.params[\"a2\"]\n", " m_break = 10 ** self.params[\"M_break\"]\n", "\n", " out = (\n", " m\n", " * (a1 * (m / 1e10) ** beta * np.exp(-((m / m_break) ** alpha)) + a2)\n", " * np.exp(-((m_1 / m) ** 0.5))\n", " )\n", " return out\n", "\n", " def _satellite_occupation(self, m):\n", " alpha = self.params[\"alpha_sat\"]\n", " beta = self.params[\"beta_sat\"]\n", " amp = 10 ** self.params[\"M_0\"]\n", " m1 = 10 ** self.params[\"M_break_sat\"]\n", " array = np.zeros_like(m)\n", " array[m >= 10**11] = 1\n", " return amp * (m / m1) ** beta * np.exp(-((m1 / m) ** alpha)) * array\n", " # return 10**8\n", "\n", " def unit_conversion(self, cosmo, z):\n", " \"A factor to convert the total occupation to a desired unit.\"\n", " A12 = 2.869e-15\n", " nu21cm = 1.42e9\n", " Const = (3.0 * A12 * const.h * const.c**3.0) / (\n", " 32.0 * np.pi * (const.m_p + const.m_e) * const.Boltzmann * nu21cm**2\n", " )\n", " Mpcoverh_3 = ((astroconst.kpc.value * 1e3) / (cosmo.H0.value / 100.0)) ** 3\n", " hubble = cosmo.H0.value * cosmo.efunc(z) * 1.0e3 / (astroconst.kpc.value * 1e3)\n", " temp_conv = Const * ((1.0 + z) ** 2 / hubble)\n", " # convert to Mpc^3, solar mass\n", " temp_conv = temp_conv / Mpcoverh_3 * astroconst.M_sun.value\n", " return temp_conv\n", "\n", " def _tracer_per_central(self, M):\n", " \"\"\"Number of tracers per central tracer source\"\"\"\n", " n_c = np.zeros_like(M)\n", " n_c[\n", " np.logical_and(\n", " 10 ** self.params[\"M_min_counts\"] <= M,\n", " 10 ** self.params[\"M_max_counts\"] >= M,\n", " )\n", " ] = 1\n", " return n_c\n", "\n", " def _tracer_per_satellite(self, M):\n", " \"\"\"Number of tracers per satellite tracer source\"\"\"\n", " n_s = np.zeros_like(M)\n", " index = np.logical_and(\n", " 10 ** self.params[\"M_min_counts\"] <= M,\n", " 10 ** self.params[\"M_max_counts\"] >= M,\n", " )\n", " n_s[index] = (M[index] / 10 ** self.params[\"M_1_counts\"]) ** self.params[\"alpha_counts\"]\n", "\n", " return n_s" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now update the halo model with this newly defined HOD:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "hm.hod_model = Spinelli19" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can check out the mean density of HI in units of $M_\\odot h^2 {\\rm Mpc}^{-3}$:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "119798016.01701438" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hm.mean_tracer_den" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And in temperature units:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "np.float64(3.7574594855492206e-05)" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hm.mean_tracer_den_unit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can easily update the parameters using:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "hm.hod_params = {\"beta\": 1.0}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And the power spectrum in length units:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAG9CAYAAAD5ixlRAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAASkNJREFUeJzt3XlYVPXiBvD3zLBvoyyyyCAuoCIKyOa+FK5laplmpmi/rNwyabnaveW9bbaamZNpZtZ1Q01N08okFTUUBUHcURGQVUQY9mVmfn9YlFclYGY4s7yf5+F54szMmZfOg/Nyzvd8v4JGo9GAiIiIyAxJxA5AREREJBYWISIiIjJbLEJERERktliEiIiIyGyxCBEREZHZYhEiIiIis8UiRERERGaLRYiIiIjMloXYAQyZWq1Gbm4uHB0dIQiC2HGIiIioCTQaDcrKyuDl5QWJpPFzPixCjcjNzYVcLhc7BhEREbVAdnY2vL29G30Oi1AjHB0dAdz+H+nk5CRyGiIiImoKpVIJuVze8DneGBahRvxxOczJyYlFiIiIyMg0ZVgLB0vfg0KhQEBAAMLDw8WOQkRERHokcPX5+1MqlZDJZCgtLeUZISIiIiPRnM9vnhEiIiIis8UiRERERGaLRYiIiIjMFosQERERmS0WISIiIjJbLEJERERktliEiIiIyGyxCBEREZHZYhEiIiIis8UiRERERGaLi66agOo6FUoq61BSVYuSyjooq+qgASARBEgEQCIRIBEEWFtI4OZoDXcnGzhY89ATERHx09AI1NSrkFFUgWtFlcgqrsC1m5XIvHn7+6LyGtTUq5u9T3srKdydbNDOyRodnO3R01uGXt4ydPNwgpUFTxQSEZF54KKrjRBj0dXaejUu5pchLacUaTklSMspxcX8MtSpGj9MUomANraWkNlZwsnGEhIBUGsAjUYDlUYDtRqorlfhhrIGZTX1992PlVSCbp6O6OUtQ//Orhjc1Q12VuzLRERkPJrz+c1PuHtQKBRQKBRQqVR6f6+aehVOZZUg4cpNJFy9iZSsEtSq7j7D42RjgY6u9ujgYo8OLnbo4GIPXxc7uDvZoI2dJRysLSAIQpPes6KmHoVlNShQVqNAWY30gnKczinF6eslKKmsw+nrpTh9vRTrj2XB2kKCQf5uGNHDA1Hd26GNnZWu/xcQERGJhmeEGqGvM0KZNyuwOzUXCVdv4uS1W3dd2mpjZ4me7WUNX4HtZfBua9vkotNSGo0G129VIfV6CZIzS7D/fAGyiisbHpdKBPTt5IInIuQY2cMDFlJeQiMiIsPTnM9vFqFG6KsI/XKuADO/PdnwvauDNfp1dkHfzi7o08kFvi52ei89TaHRaHAhvww/ncnHz2fzcSG/rOExL5kNpvb1xeQIOc8SERGRQWER0hF9FaHSqjos2n4afTvdLj+d3RwMovj8ncybFfguOQcbj2eiqLwWAGBjKcH4EG883d8Xfu6OIickIiJiEdIZMQZLG4PqOhV+OJ2HtUcycC5PCQAQBGB8cHvEDPeHd1s7kRMSEZE5YxHSERahxmk0GiRmFOOrIxnYd64AwO27zqL7dcCcoV14yYyIiETBIqQjLEJNl5pdgvd+vICEqzcBAI42Fpg9pAtm9PeFjaVU5HRERGROWIR0hEWoeTQaDQ5duoH3frzQMLDau60t3hoXiKFd24mcjoiIzAWLkI6wCLWMSq3B9yk5+Ojni8gtrQYAPNTLE4sfDkA7JxuR0xERkalrzuc3J4IhnZNKBDza2xu/xAzGzIEdIZUI2HM6Dw9+fAj/PZYJtZrdm4iIDAOLEOmNvbUF/vlQAL6f0x9B3jKU1dTj9Z1n8NgXv+HiX+YkIiIiEguLEOldYHsZts/uj3+PCYCDtQVOZZVgzGdH8PnBy6i/x3IiRERErYVFiFqFVCJgev+O+CVmEB7s1g61KjU++OkiJnyRgMuF5WLHIyIiM8UiRK3KU2aLNdFh+OjxIDjaWCAluwSjlx/Gl/FXoeLYISIiamUsQtTqBEHAhFBv7FswCIP83VBbr8Y7e89j0qoEZBRViB2PiIjMCIsQicZTZotvZoRjyaM9YW8lxcnMWxj1aTzWHc3gnWVERNQqWIRIVIIgYHKED35eMAj9Oruguk6Nf+8+hylrjiO7uFLseEREZOJYhMggeLe1w/r/i8SbY3vA1lKKhKs3MXJZPDYezwLn/CQiIn1hESKDIZEImNbXFz/OH4hw37aoqFXhtR1pmLY2kWeHiIhIL8yiCFVWVqJDhw54+eWXxY5CTeDrao/Nz/bFvx7qDmsLCQ6nF2HYJ4ewOv4K5x0iIiKdMosi9M4776BPnz5ix6BmkEoEPDOwE356cRD6dro9dujdvRcwVnEUaddLxY5HREQmwuSLUHp6Oi5cuIBRo0aJHYVaoKOrPTbOjMQHE3pBZmuJs7lKjFUcwVs/nENFTb3Y8YiIyMgZdBGKj4/HmDFj4OXlBUEQsHPnzrueo1Ao4OvrCxsbG0RGRiIxMfGOx19++WUsWbKklRKTPgiCgIlhcsS9NBiPBHlBrQG+OpKBoR8dxNaT2bzVnoiIWsygi1BFRQWCgoKgUCju+XhsbCxiYmKwePFiJCcnIygoCCNGjEBhYSEA4Pvvv4e/vz/8/f1bMzbpiauDNZZPDsHXM8Ihd7ZFYVkNXtl2GmNWHEHClZtixyMiIiMkaIzk3mRBELBjxw6MGzeuYVtkZCTCw8OxYsUKAIBarYZcLse8efOwcOFCLFq0COvXr4dUKkV5eTnq6urw0ksv4Y033rjne9TU1KCmpqbhe6VSCblcjtLSUjg5Oen156PmqalX4ZvfruGzuMso+/0S2fAAdywa3R0dXe1FTkdERGJSKpWQyWRN+vw26DNCjamtrUVSUhKioqIatkkkEkRFRSEhIQEAsGTJEmRnZ+PatWv46KOPMHPmzPuWoD+eL5PJGr7kcrnefw5qGWsLKZ4d1BkHXxmCqX06QCoRsO9cAYYtPYR/7khDbkmV2BGJiMgIGG0RKioqgkqlgru7+x3b3d3dkZ+f36J9Llq0CKWlpQ1f2dnZuohKeuTiYI23xgXip/kDMaSrG+rVGmw4noUhHx7Ev3edRaGyWuyIRERkwCzEDtBapk+f/rfPsba2hrW1tf7DkM75uTti3YwIHLt6E0t/uYTEjGKs++0aNiVmYWqfDnh+SGe4OvDYEhHRnYz2jJCrqyukUikKCgru2F5QUAAPDw+RUpHY+nRyQeyzfbDhmUj09mmDmno11hzJwKAPDuDT/emorOUt90RE9CejLUJWVlYIDQ1FXFxcwza1Wo24uDj07dtXq30rFAoEBAQgPDxc25gkAkEQ0L+LK76b1Q/fPB2BIG8ZKmtV+GT/JQz58CC2nMiGirfcExERDPyusfLycly+fBkAEBISgqVLl2Lo0KFwdnaGj48PYmNjER0djVWrViEiIgLLli3Dli1bcOHChbvGDrVEc0adk+HSaDTYk5aH93+6gOzi24Oou3k44rXR3THI303kdEREpGvN+fw26CJ08OBBDB069K7t0dHRWLduHQBgxYoV+PDDD5Gfn4/g4GAsX74ckZGROnl/FiHTUlOvwre/ZeKzX9OhrL59iWxED3d8MCEIMltLkdMREZGumEwREhuLkGm6VVGLFQcu49uEa6hTadDBxQ5fPBWK7p48xkREpsAs5hHSJ44RMm1t7a3w+sMB+G5WP7RvY4vMm5UY//lR7Dh1XexoRETUynhGqBE8I2T6iitqMX/zKRxOLwIATO3TAa8/HAArC/6NQERkrHhGiKiJnO2tsG5GBF54oAsA4L/HMjFpdQInYiQiMhMsQmT2pBIBMcO7Yu30MDjZWOBUVgmmrU1EWXWd2NGIiEjPWISIfvdAN3fsmjsAbo7WuJBfhtkbklGnUosdi4iI9IhF6B44WNp8+bra46voMNhaSnE4vQj/2nEGHEZHRGS6OFi6ERwsbb7izhdg5rcnodYALw/3x9wH/MSORERETcTB0kRaerC7O/7zSA8AwEf7LmHnqRyRExERkT6wCBHdx9S+vnh2UCcAwCvbUpFw5abIiYiISNdYhIgasXBkNzzU0xN1Kg2e++9JXCuqEDsSERHpEIvQPXCwNP1BIhHw8cQg9PZpA2V1PeZvPsU7yYiITAgHSzeCg6XpD7klVRi5LB7K6nrMHtIZr47sJnYkIiK6Dw6WJtIxrza2eO+xXgCAlYeu4LcrRSInIiIiXWARImqi0T098US4HBoNEBObilsVtWJHIiIiLbEIETXDG2MC0MnVHvnKaizcfpqTLRIRGTkWIaJmsLOywPLJIbCUCvj5bAE2JWaLHYmIiLTAIkTUTIHtZXh1xO3B0m/+cBaXC8tETkRERC3FInQPvH2e/s7/DeiIgX6uqK5TY+7GU6ioqRc7EhERtQBvn28Eb5+nxhSWVWP0p0dQVF6DB7q1w5fTwiCVCGLHIiIye7x9nqgVtHO0wZfTQmFtIcGvFwrx1g/nxI5ERETNxCJEpIUQn7ZYNikYALDut2v4+miGuIGIiKhZWISItDSqpycWjro9ePqtH85h/7kCkRMREVFTsQgR6cBzgzphcoQcag0wb9MpnMkpFTsSERE1AYsQkQ4IgoA3xwZioJ8rqupUeHrdCeSWVIkdi4iI/gaLEJGOWEolUEzpDX93BxSW1WCs4ih+u8w1yYiIDBmL0D1wHiFqKScbS6ydHg6/dg64UVaDKV8dx8f7LqJepRY7GhER3QPnEWoE5xGilqqqVeE/u89i84nbS3BE+Drj08nB8JTZipyMiMj0cR4hIpHZWknx3mO9sHxyCBysLZB4rRijPj3MO8qIiAwMixCRHj0S5IUf5g1Az/YylFTW4ZlvT2J1/BWxYxER0e9YhIj0zNfVHt/N6ofp/XwBAO/uvYBl+y+BV6WJiMTHIkTUCqwsJPj3Iz3wyoiuAIBl+9Px3o8XWIaIiETGIkTUiuYM7YI3Hg4AAKyKv4o3vj8LtZpliIhILCxCRK3s6QEdseTRnhAE4L/HMvHqd6ehYhkiIhIFixCRCCZH+GDpxCBIJQK2JV3Hi7EpvExGRCQCFiEikYwP8YbiyRBYSgXsTs3FlpPZYkciIjI7LEL3wJmlqbWMDPTEqyNur1z/zp7zKCyrFjkREZF54czSjeDM0tQa6lVqjPv8KM7kKPFwL0+seLK32JGIiIwaZ5YmMiIWUgnee7QXpBIBP5zOw68XOPs0EVFrYREiMgCB7WX4vwEdAQD/2nEG5TX1IiciIjIPLEJEBuLFKD/InW2RW1qNj/ddFDsOEZFZYBEiMhB2VhZ4Z1xPAMC6364hJbtE3EBERGaARYjIgAzyd8P4kPbQaICF351GnUotdiQiIpPGIkRkYP71UHe0tbPEhfwyrI6/KnYcIiKTxiJEZGBcHKzx+u/rkS395RISM4pFTkREZLpYhIgM0PiQ9hgb7AWVWoPZG5JRoOREi0RE+sAiRGSABEHAkkd7oqu7I4rKazB7QzJq6zleiIhI11iEiAyUnZUFvpgaCkdrCyRl3sK7e8+LHYmIyOSwCBEZsI6u9lg6KRjA7Vvqd57KETcQEZGJYREiMnDDAtwxd2gXAMDC7adxPk8pciIiItPBIkRkBBYM88dAP1dU16nx/PoklFbWiR2JiMgksAjdg0KhQEBAAMLDw8WOQgQAkEoELH8iBO3b2CLzZiWiv06EsppliIhIW4JGo9GIHcJQKZVKyGQylJaWwsnJSew4RDibW4onvzyO0qo69PKW4b9PR0JmZyl2LCIig9Kcz2+eESIyIj28ZNg4MxJt7Sxx+nopnlxzDLcqasWORURktFiEiIxMDy8ZNj3bBy72Vjibq8STa47jZnmN2LGIiIwSixCREerm4YTNz/aBq4M1zucp8eSXx1HEMkRE1GwsQkRGys/dEZuf7YN2jta4WFCGyauPcQA1EVEzsQgRGbEu7RwQ+1xfeDjZIL2wHB/+dFHsSERERoVFiMjIdXS1x9KJQQCA9cczkZR5S+RERETGg0WIyAT06+KKCaHe0GiA17ancYFWIqImYhEiMhH/HN0dzvZWuFhQhi8PXxU7DhGRUWARIjIRbe2t8PrD3QEAn8al41pRhciJiIgMH4sQkQkZF9weA/1cUVuvxj93poETxxMRNY5FiMiECIKAt8cFwtpCgqOXb2J7co7YkYiIDBqLEJGJ6eBij/lRfgCAt/ecQzGX4CAiui8WISITNHNgJ3TzcMStyjq8ufssL5EREd0HixCRCbKUSrDk0Z4QBGBnSi6+TcgUOxIRkUFiESIyUSE+bbFwZDcAwJs/nMPh9BsiJyIiMjwsQkQm7NlBnfBo7/ZQqTWYsyEZV2+Uix2JiMigsAgRmTBBEPDu+J7o7dMGyup6PPPNSZRWcmFWIqI/mHQRKikpQVhYGIKDgxEYGIgvv/xS7EhErc7GUoovpobCS2aDq0UVmLspGfUqLsFBRAQAgsaEbydRqVSoqamBnZ0dKioqEBgYiJMnT8LFxaVJr1cqlZDJZCgtLYWTk5Oe0xLp19ncUkxYmYCqOhVm9PfF4jE9xI5ERKQXzfn8NukzQlKpFHZ2dgCAmpoaaDQa3kZMZquHlwyfTLq9Sv3XR69hy8lskRMREYnPoItQfHw8xowZAy8vLwiCgJ07d971HIVCAV9fX9jY2CAyMhKJiYl3PF5SUoKgoCB4e3vjlVdegaurayulJzI8IwM9ETPMHwDw711nkXWzUuRERETiMugiVFFRgaCgICgUins+Hhsbi5iYGCxevBjJyckICgrCiBEjUFhY2PCcNm3aIDU1FRkZGdi4cSMKCgpaKz6RQZo7tAsiOjqjslaFl7elQq3mWVIiMl8GXYRGjRqFt99+G+PHj7/n40uXLsXMmTMxY8YMBAQE4IsvvoCdnR3Wrl1713Pd3d0RFBSEw4cP3/f9ampqoFQq7/giMjUSiYCPHw+CvZUUiRnFWHs0Q+xIRESiMegi1Jja2lokJSUhKiqqYZtEIkFUVBQSEhIAAAUFBSgrKwMAlJaWIj4+Hl27dr3vPpcsWQKZTNbwJZfL9ftDEIlE7myHfz0cAAD44OeLuFxYJnIiIiJxGG0RKioqgkqlgru7+x3b3d3dkZ+fDwDIzMzEwIEDERQUhIEDB2LevHno2bPnffe5aNEilJaWNnxlZ3MwKZmuJ8LlGNLVDbX1asRsSUUdb6knIjNkIXYAfYqIiEBKSkqTn29tbQ1ra2v9BSIyIIIg4P3HemH4J/E4fb0UKw9ewQsP+okdi4ioVRntGSFXV1dIpdK7Bj8XFBTAw8NDpFRExsXdyQZvjr09n9DyuHScySkVORERUesy2iJkZWWF0NBQxMXFNWxTq9WIi4tD3759tdq3QqFAQEAAwsPDtY1JZPAeCfLC6J4eqFdrELMlBdV1KrEjERG1GoMuQuXl5UhJSWm4vJWRkYGUlBRkZWUBAGJiYvDll1/im2++wfnz5zFr1ixUVFRgxowZWr3vnDlzcO7cOZw4cULbH4HI4AmCgLfH9YSrgxUuFZTjvR8viB2JiKjVGPQYoZMnT2Lo0KEN38fExAAAoqOjsW7dOkyaNAk3btzAG2+8gfz8fAQHB+Onn366awA1ETXO2d4KHz4ehBlfn8C6366hTycXjAzkJWYiMn0mvdaYtrjWGJmbJXvPY1X8VTjZWGDPCwMhd7YTOxIRUbNxrTEtcYwQmauXR3RFiE8bKKvrMW/TKd5ST0Qmj2eEGsEzQmSOsosr8dDyw1BW1+O5QZ2waHR3sSMRETVLcz6/WzRGaNeuXc1+zbBhw2Bra9uStyOiViR3tsMHE4Lw/PokrIq/ij6dXDC0WzuxYxER6UWLzghJJM27oiYIAtLT09GpU6fmvpWoeEaIzNm/d53Fut+uoa2dJfbOHwhPGf+QISLj0CpjhPLz86FWq5v0ZWfHAZdExmbR6G4IbO+EW5V1mL85BSquUk9EJqhFRSg6OrpZl7meeuopozqjwsHSRIC1hRQrJveGg7UFEjOK8TVXqSciE8TB0o3gpTEiYFNiFhZtT4O1hQR75w9EZzcHsSMRETVKr5fGbt26heLiYgDAjRs3sH37dpw9e7ZlSYnI4D0RLsdAP1fU1KvxytZUXiIjIpPSrCK0Zs0ahIaGIiwsDCtXrsT48eMRFxeHJ554AmvWrNFXRiIS0R+r1DtaWyA5qwRfHbkqdiQiIp1p1qWxXr164fjx46iqqoKPjw8yMjLg5uaG0tJSDB48uGFNMFPBS2NEf9pyIhuvfncaVhYS7H1hALq0cxQ7EhHRPent0piFhQVsbW3h7OyMLl26wM3NDQAgk8kgCELLExsYDpYmutvjYd4Y0tUNtfVqvLT1NOo56zQRmYBmFSGpVIrq6moAwKFDhxq2l5eX6zaVyLj6PNHdBEHAkkd7wtHGAqnZJfjyMO8iIyLj16witH//flhbWwO4fRboD5WVlVi9erVukxGRwfGU2eKNhwMAAJ/8cgnpBWUiJyIi0k6zitD/XgLLz88HALRr146XkYjMxIRQbzzQrR1qVWrEbEnlwqxEZNS0Wn1++PDhuspBREbij0tkMltLpOWUYnlcutiRiIhaTKsixLkYicyTu5MN3h3fEwCgOHAZSZnFIiciImoZrYqQKd0pRkTN81AvTzwa0h5qDbAgNhXlNfViRyIiajatipCp4u3zRE3z77E90L6NLbKKK/HW7nNixyEiajYWoXvg7fNETeNkY4mPJwZBEIDYk9n4+Wy+2JGIiJpFqyIklUp1lYOIjFSfTi54dmAnAMCi7WkoLKsWORERUdNpVYROnTqlqxxEZMRihvuju6cTiitq8Y9tp3kjBREZDV4aIyKtWVtIsWxSMKwsJDhw8QY2JmaJHYmIqEksdLGTuLg4xMXFobCwEGr1nZOrrV27VhdvQUQGrquHI14d0RVv7zmPd/acx8AubvBxsRM7FhFRo7Q+I/Sf//wHw4cPR1xcHIqKinDr1q07vojIfDzdvyMiOzqjslaFl7emQqXmJTIiMmyCRsuL+Z6envjggw8wdepUXWUyGEqlEjKZDKWlpXBychI7DpFRyC6uxMhl8aioVeGfo7tj5qBOYkciIjPTnM9vrc8I1dbWol+/ftruxqBwHiGilpM72+Ffvy/M+uG+i1yYlYgMmtZF6JlnnsHGjRt1kcVgcB4hIu08ES7HkK5uqK3nwqxEZNhaNFg6Jiam4b/VajVWr16N/fv3o1evXrC0tLzjuUuXLtUuIREZHUEQ8P5jvTD8k3ik5ZRCceAyXozyFzsWEdFdWlSE/nf+oODgYADAmTNn7tjOtciIzJe7kw3eHNsD8zenYMWvl/FgN3f09JaJHYuI6A5aD5Y2ZRwsTaQdjUaDuRtPYU9aHvzaOWD3vAGwseSM9ESkX606WJqI6H4EQcBb4wLh6mCN9MJyfLL/ktiRiIju0OwidOvWLRQXFwMAbty4ge3bt+Ps2bM6D0ZEpsHZ3gpLHu0JAPgy/iqSMjm/GBEZjmYVoTVr1iA0NBRhYWFYuXIlxo8fj7i4ODzxxBNYs2aNvjISkZEbFuCOR3u3h1oDvLw1FVW1KrEjEREBaOYYoV69euH48eOoqqqCj48PMjIy4ObmhtLSUgwePBgpKSl6jNr6OEaISHdKK+swfNkhFChrMKO/LxaP6SF2JCIyUXobI2RhYQFbW1s4OzujS5cucHNzAwDIZDLeIUZEjZLZWeL9x3oBAL4+eg3Hrt4UORERUTOLkFQqRXV1NQDg0KFDDdvLy8t1m0pknFmaSD+GdG2HJ8LlAIBXtqWioqZe5EREZO6adWnsj1NM/3v2p7CwEJmZmSZXHHhpjEj3yqrrMHLZYeSUVGFKpA/eGd9T7EhEZGL0dmnsfpfA2rVrZ3IliIj0w9HGEh9OuH2JbMPxLBxOvyFyIiIyZ1rPI7RkyRKsXbv2ru1r167F+++/r+3uicgE9eviiml9OwAAXt12GqVVdSInIiJzpXURWrVqFbp163bX9h49euCLL77QdvdEZKIWjuqGDi52yCutxls/nBM7DhGZKa2LUH5+Pjw9Pe/a7ubmhry8PG13T0Qmys7KAh89HgRBALYlXccv5wrEjkREZkjrIiSXy3H06NG7th89ehReXl7a7p6ITFi4rzNmDuwEAFi0PQ3FFbUiJyIic6N1EZo5cyZefPFFfP3118jMzERmZibWrl2LBQsWYObMmbrISEQmLGaYP/zaOaCovAavf39G7DhEZGYstN3BK6+8gps3b2L27Nmora2FRqOBra0t/vGPf2DRokW6yEhEJszGUoqlE4Mx7vOj2HM6DyN75GJMEM8mE1HraNY8Qo0pLy/HuXPnYGtrC39/f1hbW+tit6LiPEJEreeTXy7h07h0tLGzxL4XB6Gdk43YkYjISOltHqH7+eqrr9CnTx8MHDgQYWFhCA0N5SKsRNQscx/ogh5eTiiprMOi7WnQ0d9oRESN0roIvfHGG5g/fz7GjBmDrVu3YuvWrRgzZgwWLFiAN954QxcZicgMWEolWDoxGFZSCeIuFGLryetiRyIiM6D1pTE3NzcsX74ckydPvmP7pk2bMG/ePBQVFWkVUEy8NEbU+r44dAXv/XgBDtYW+HH+QMid7cSORERGplUvjdXV1SEsLOyu7aGhoaiv54KKRNQ8Mwd2QliHtiivqcfLW1OhVvMSGRHpj9ZFaOrUqVi5cuVd21evXo0pU6Zou3siMjNSiYCPJwbBzkqK4xnFWHs0Q+xIRGTCtL59Hrg9WHrfvn3o06cPAOD48ePIysrCtGnTEBMT0/C8pUuX6uLt9E6hUEChUEClUokdhcgsdXCxx78eCsBrO9Lwwc8XMdjfDX7ujmLHIiITpPUYoaFDhzbtjQQBv/76qzZv1eo4RohIPBqNBjPWncDBizcQ2N4JO2b3h6VUJze6EpGJa87nt87mETJFLEJE4ipUVmP4sniUVNbhhQf9EDPMX+xIRGQEWn0eISIifWjnZIO3xwUCABQHLiMlu0TcQERkclo8Rujpp59u0vPWrl3b0rcgIsLDvbyw72wBdqXmIiY2BXteGAhbK6nYsYjIRLS4CK1btw4dOnRASEgIZ4AlIr16c2wPHM+4iatFFXjvx/P4z9hAsSMRkYlocRGaNWsWNm3ahIyMDMyYMQNPPfUUnJ2ddZmNiAgA0MbOCh9OCMK0tYn4JiETD3R3x2B/N7FjEZEJaPEYIYVCgby8PLz66qvYvXs35HI5Jk6ciJ9//plniIhI5wb5u2F6P18AwCtbU3GrolbcQERkErQaLG1tbY3Jkyfjl19+wblz59CjRw/Mnj0bvr6+KC8v11VGIiIAwD9GdkNnN3sUltXgtR1cmJWItKezu8YkEgkEQYBGo+FEhESkF7ZWUnz6RAgsJAJ+PJOP7ck5YkciIiOnVRGqqanBpk2bMGzYMPj7+yMtLQ0rVqxAVlYWHBwcdJWRiKhBYHsZFvw+n9DiXWeRXVwpciIiMmYtLkKzZ8+Gp6cn3nvvPTz88MPIzs7G1q1bMXr0aEgknJ6IiPTn+cGdGxZmfWlLKlRcmJWIWqjFM0tLJBL4+PggJCQEgiDc93nbt29vcTixcWZpIsOVXVyJkcviUVGrwj9GdsOsIZ3FjkREBqI5n98tvn1+2rRpjRYgIiJ9kjvbYfEjPfDqttNY+stFDPRzRWB7mdixiMjIcK2xRvCMEJFh02g0mL0hGT+eyUcnN3v8MG8A7Kxa/PcdEZkIrjVGRGZBEAQsebQnPJxscPVGBd764bzYkYjIyLSoCJ0+fRpqtbrJzz979izq6+tb8lZERI1qY2eFpRODIAjApsQs/Hw2X+xIRGREWlSEQkJCcPPmzSY/v2/fvsjKymrJWxER/a1+XVzx7KBOAICF351GgbJa5EREZCxadDFdo9Hg9ddfh52dXZOeX1vLqfCJSL9eGtYVRy8X4UyOEi9tScW3T0dAIuENHUTUuBYVoUGDBuHixYtNfn7fvn1ha2vbkrfSSnZ2NqZOnYrCwkJYWFjg9ddfx+OPP97qOYhI/6wsJFg2KQQPf3YYRy4X4asjGZj5+1kiIqL7Mem7xvLy8lBQUIDg4GDk5+cjNDQUly5dgr29fZNez7vGiIzPxuNZeG1HGiylAnbM7s9b6onMEO8a+52npyeCg4MBAB4eHnB1dUVxcbG4oYhIryZHyDEswB11Kg1e2HwKlbW8UYOI7s+gi1B8fDzGjBkDLy8vCIKAnTt33vUchUIBX19f2NjYIDIyEomJiffcV1JSElQqFeRyuZ5TE5GYBEHA+4/1gruTNa7eqMC/d50VOxIRGTCDLkIVFRUICgqCQqG45+OxsbGIiYnB4sWLkZycjKCgIIwYMQKFhYV3PK+4uBjTpk3D6tWrWyM2EYnM2d4KyyaFQBCALSevY1dqrtiRiMhAGc0YIUEQsGPHDowbN65hW2RkJMLDw7FixQoAgFqthlwux7x587Bw4UIAQE1NDYYNG4aZM2di6tSpjb5HTU0NampqGr5XKpWQy+UcI0RkpD7edxGf/XoZDtYW2PvCQPi4NO1OVyIybmYxRqi2thZJSUmIiopq2CaRSBAVFYWEhAQAt2/znz59Oh544IG/LUEAsGTJEshksoYvXkYjMm7zH/RrWKV+3uZTqFM1fSJYIjIPOilCdXV1yM7OxsWLF1ttMHJRURFUKhXc3d3v2O7u7o78/Nszyx49ehSxsbHYuXMngoODERwcjLS0tPvuc9GiRSgtLW34ys7O1uvPQET6ZSGVYNkTwXCysUBqdgk+2tf0aT+IyDy0eHXCsrIyrF+/Hps3b0ZiYiJqa2uh0WggCAK8vb0xfPhwPPvsswgPD9dl3mYZMGBAs5YCsba2hrW1tR4TEVFr825rhw8m9MLz65Ox6tBV9O/sikH+bmLHIiID0aIzQkuXLoWvry++/vprREVFYefOnUhJScGlS5eQkJCAxYsXo76+HsOHD8fIkSORnp6u69xwdXWFVCpFQUHBHdsLCgrg4eGh8/cjIuM1MtATUyJ9AAAxW1Jxo6zmb15BROaiRWeETpw4gfj4ePTo0eOej0dERODpp5/GypUrsW7dOhw+fBh+fn5aBf1fVlZWCA0NRVxcXMMAarVajbi4OMydO1erfSsUCigUCqhUKh0kJSJD8PrDATh57RYuFpRhQWwKvnk6AlIuwUFk9gz6rrHy8nJcvnwZwO2FXpcuXYqhQ4fC2dkZPj4+iI2NRXR0NFatWoWIiAgsW7YMW7ZswYULF+4aO9QSnFmayLSkF5ThkRVHUVWnwkvD/DHvQd3+gUZEhqFV7xrbtGnTfR975ZVXtNr3yZMnERISgpCQEABATEwMQkJC8MYbbwAAJk2ahI8++ghvvPEGgoODkZKSgp9++kknJYiITI+fuyPeHhcIAPhk/yUkXLkpciIiEpvWZ4TatGmDTZs2YdSoUXdsX7BgATZv3oy8vDytAoqJZ4SITNMrW1OxNek63BytsfeFgXBz5E0SRKakVc8IbdiwAZMnT8aRI0cats2bNw9btmzBgQMHtN29KBQKBQICAkS9442I9OfNsYHwd3fAjbIaLIhNgUptsCMEiEjPdDJGaOPGjZg7dy5++eUXfPXVV/j+++9x4MAB+Pv76yKjaHhGiMh0cbwQkelqzud3i+cR+qsnn3wSJSUl6N+/P9zc3HDo0CF06dJFF7smItKLP8YLvbQ1FZ/sv4QwX2f07ewidiwiamUtKkIxMTH33O7m5obevXvj888/b9i2dOnSliUjItKzx0K9cezqTWxNuo4XNp/CnhcGoJ2jjdixiKgVtagInTp16p7bu3TpAqVS2fC4IHCODiIybG+ODUTq9RJcKijHvI2nsOGZSFhIjXYZRiJqJoOeR0gsf51Q8dKlSxwjRGTirtwox9gVR1FeU4/nBnXCotHdxY5ERFrQ+11jWVlZzXp+Tk5OS95GNHPmzMG5c+dw4sQJsaMQUSvo7OaADyf0AgCsir+Kn84Y77QfRNQ8LSpC4eHheO655xotCqWlpfjyyy8RGBiI7777rsUBiYhaw6ienpg5sCMA4OWtp3H1RrnIiYioNbRojNC5c+fwzjvvYNiwYbCxsUFoaCi8vLxgY2ODW7du4dy5czh79ix69+6NDz74AKNHj9Z1biIinfvHyG5IvV6KxIxizFqfjB1z+sHOSic31xKRgdJqjFBVVRX27NmDI0eOIDMzE1VVVXB1dUVISAhGjBiBwMBAXWZtdZxHiMj8FCqr8dBnR3CjrAbjgr3wyaRg3vhBZGSa8/nNwdKNYBEiMk+JGcWY/OUxqNQavDm2B6b19RU7EhE1Q6sssaFSqfD++++jf//+CA8Px8KFC1FVVdXS3RkULrFBZN4iOjpj0ahuAIA3d59DYkaxyImISF9aXITeffddvPbaa3BwcED79u3x6aefYs6cObrMJhreNUZE/zegIx7u5Yl6tQazNyQht8Q0/tAjoju1uAh9++23+Pzzz/Hzzz9j586d2L17NzZs2AC1Wq3LfEREohAEAR9M6IXunk4oKq/Fc/9NQnWdSuxYRKRjLS5CWVlZd9wNFhUVBUEQkJubq5NgRERis7OywOqpoXC2t0JaTikWbU8Dh1USmZYWF6H6+nrY2Ny5Jo+lpSXq6uq0DkVEZCjkznZY8WQIpBIBO07l4KsjGWJHIiIdavEEGRqNBtOnT4e1tXXDturqajz//POwt7dv2LZ9+3btEhIRiaxfZ1f866Hu+M/uc3h373l083DCAD9XsWMRkQ60uAhFR0ffte2pp57SKoyh+OtaY0READC9ny/O5iqxLek65mxMxq65/dHBxf7vX0hEBo3zCDWC8wgR0V9V16kwafUxpGaXoEs7B3w3qx9ktpZixyKi/9Eq8wgREZkbG0spVk8NhYeTDS4XlmPuxmTUqXinLJExYxEiImoGdycbrIkOg62lFIfTi/DvXWd5JxmREWMRIiJqpsD2MiyfHAJBADYcz8LXR6+JHYmIWohFiIioBYYFuOO1Ud0BAG/vOYdfLxSInIiIWoJFiIiohZ4Z2BGTI+RQa4B5G0/hfJ5S7EhE1EwsQkRELSQIAt4cG4h+nV1QUavC/607gQJltdixiKgZWITugavPE1FTWUolWDklFJ3c7JFbWo3otYlQVnOGfSJjwXmEGsF5hIioqbKLK/Hoyt9wo6wGfTo5Y92MCNhYSsWORWSWOI8QEVErkzvbYd2McDhYW+DY1WLEbEmBSs2/M4kMHYsQEZGO9PCSYfXUUFhKBexNy8ebuznHEJGhYxEiItKhfl1csXRiMADgm4RMrDx0RdxARNQoFiEiIh0bE+SFNx4OAAB88NNFbD2ZLXIiIrofFiEiIj14ekBHPDe4EwDgH9+dxt60PJETEdG9sAgREenJP0Z0w+Oh3lBrgBc2nULcec4+TWRoWISIiPREIhHw3mO98EiQF+rVGsxan4zD6TfEjkVEf8EiRESkR1KJgI8nBmFED3fUqtSY+e1JHL96U+xYRPQ7FiEiIj2zlEqwfHIIhnR1Q3WdGk+vO4HkrFtixyIisAjdE5fYICJds7aQ4ounQhvWJYtem4gzOaVixyIye1xioxFcYoOIdK2yth7TvkrEycxbcLKxwDdPRyDEp63YsYhMCpfYICIyUHZWFvh6RjjCOrSFsroeT605zjFDRCJiESIiamWONpb49v8i/rxM9nUi4i/xbjIiMbAIERGJwM7KAmunh2Po7wOon/nmJH45x3mGiFobixARkUhsLKVYNTUMI3t4oFalxqz1SfjhdK7YsYjMCosQEZGIrCwkWPFkCMYF35508YVNp7DxeJbYsYjMBosQEZHILKQSfDwxGJMj5FBrgNd2pOHDny+AN/US6R+LEBGRAZBKBLw7vifmP+gHAFAcuIKYLamorVeLnIzItLEIEREZCEEQsGCYP95/rCekEgE7TuVgxrpEKKvrxI5GZLJYhIiIDMykcB98FR0Geyspjl6+iYlfJCCvtErsWEQmiUWIiMgADenaDrHP9YWbozUu5JdhnOIoUrJLxI5FZHJYhIiIDFRgexm2z+qHLu0cUKCswcRVCdhyMlvsWEQmhUWIiMiAyZ3tsGN2P0R1d0dtvRqvbjuNxd+fQZ2Kg6iJdIFF6B64+jwRGRJHG0usnhqKF6Nu31H2TUImpqw5jqLyGpGTERk/rj7fCK4+T0SGZt/ZfMRsSUV5TT28ZDb4Ymooenm3ETsWkUHh6vNERCZqeA8P7JzTD51c7ZFbWo3HVv6GNYevQq3m37RELcEiRERkZLq0c8TOuf0xooc76lQavL3nPGasO8FLZUQtwCJERGSEnGws8cVToXh7XCCsLSQ4dOkGRi47jMPpN8SORmRUWISIiIyUIAh4qk8H7Jo7AH7tHFBUXoOpXyViyY/nuTQHUROxCBERGbmuHo7YNXcAnoz0AQCsOnQVj6w4gjM5pSInIzJ8LEJERCbA1kqKd8f3xMopvdHWzhIX8sswVnEUH/18ETX1KrHjERksFiEiIhMyqqcn9i0YjNE9PaBSa7DiwGWM+ewIUrk8B9E9sQgREZkYN0drfD4lFJ9P6Q0XeytcKijH+M+PYsmP51FZWy92PCKDwiJERGSiRvf0xC8xg/FIkBfUmttjh4YtjcdPZ/LAuXSJbmMRIiIyYc72Vlg+OQRrpoWhfRtb5JRU4fn1yYj++gSu3igXOx6R6FiEiIjMQFSAO/bHDMa8B7rASipB/O/zDn348wVeLiOzxrXGGsG1xojIFGUUVeA/u8/i4MXbky96ONkgZrg/HuvtDalEEDkdkfaa8/nNItQIFiEiMlUajQa/nCvAf3afQ05JFQCgq7sjFo7uhiH+bhAEFiIyXixCOsIiRESmrrpOhf8mZGLFgcsoraoDAPTr7IJFo7qjp7dM5HRELcMipCMsQkRkLkor6/D5wcv4+rdrDctzjO7pgfkP+qOrh6PI6Yiah0VIR1iEiMjcXL9ViaX7LmFHSg40GkAQbt+GP/9BP/i7sxCRcWAR0hEWISIyVxfzy7A8Lh170vIA3C5ED/1eiPxYiMjANefz2+Rvnx8/fjzatm2LCRMmiB2FiMhodPVwhGJKb/z04kCM7ukBjQb44XQehi+Lx7PfnkRy1i2xIxLphMmfETp48CDKysrwzTffYNu2bc16Lc8IERHddj5PiU/3p+Ons/kN2yI7OuP5IZ15lxkZHJ4R+oshQ4bA0ZGncYmItNHd0wlfTA3F/phBeDzUG5ZSAcczijHj6xMY9elh7Dh1vWGQNZExMegiFB8fjzFjxsDLywuCIGDnzp13PUehUMDX1xc2NjaIjIxEYmJi6wclIjITXdo54sPHgxD/6lDMHNgR9lZSXMgvw4LYVAx4/1es+DUdN8trxI5J1GQGXYQqKioQFBQEhUJxz8djY2MRExODxYsXIzk5GUFBQRgxYgQKCwtbOSkRkXnxlNninw8F4LeFD+KVEV3RztEahWU1+GjfJfR971f8Y9tpXMhXih2T6G8ZzRghQRCwY8cOjBs3rmFbZGQkwsPDsWLFCgCAWq2GXC7HvHnzsHDhwobnHTx4ECtWrPjbMUI1NTWoqfnzLxmlUgm5XM4xQkREf6O2Xo29aXlYezQDp6+XNmzv08kZ0/r6YliAOyylBv23N5kQsxgjVFtbi6SkJERFRTVsk0gkiIqKQkJCQov2uWTJEshksoYvuVyuq7hERCbNykKCcSHt8f2c/tj2fF+M7ukBiQAcu1qM2RuSMeD9X7Fs/yUUKKvFjkp0B6MtQkVFRVCpVHB3d79ju7u7O/Lz/7yrISoqCo8//jj27t0Lb2/vRkvSokWLUFpa2vCVnZ2tt/xERKZIEASE+Trj8ymhOPKPBzB3aBe4OlihQFmDZfvT0e+9XzF7QxKOXi6CWm0UFyTIxFmIHUDf9u/f3+TnWltbw9raWo9piIjMh1cbW7w8oiteeNAPP53Nx/qETCReK8betHzsTctHBxc7PBHugwmh3nBz5L+9JA6jLUKurq6QSqUoKCi4Y3tBQQE8PDxESkVERP/LykKCR4K88EiQFy7kK7HhWBZ2nspB5s1KvP/TBSz95SKGB3hgcoQP+nV2gUTCOYmo9RjtpTErKyuEhoYiLi6uYZtarUZcXBz69u2r1b4VCgUCAgIQHh6ubUwiIvqLbh5OeGtcII7/80F8MKEXguVtUKfSYE9aHp766jgGfXgAy+PSkVdaJXZUMhMGfddYeXk5Ll++DAAICQnB0qVLMXToUDg7O8PHxwexsbGIjo7GqlWrEBERgWXLlmHLli24cOHCXWOHWoIzSxMR6d+5XCU2JWZhZ0oOyqrrAQASARjs74ZJ4XI82J13nFHzmMyiqwcPHsTQoUPv2h4dHY1169YBAFasWIEPP/wQ+fn5CA4OxvLlyxEZGamT92cRIiJqPVW1Kvx0Ng+bE7NxPKO4YburgxXGh7THxDA5F3ylJjGZIiQ2FiEiInFkFFVgy8lsbEu6jhtlf87vFuLTBhPD5Hi4lyccbSxFTEiGjEVISwqFAgqFAiqVCpcuXWIRIiISSb1KjYMXb2DLyWz8eqEQ9b/fcm9rKcWonh6YGCZHZEdnLvpKd2AR0hGeESIiMhw3ymqw49R1xJ7IxpUbFQ3bO7jY4fFQbzwW6g1Pma2ICclQsAjpCIsQEZHh0Wg0SM4qwbakbOxOzUN5zZ8DrAf6uWFimBxRAe1gbSEVOSmJhUVIR1iEiIgMW2VtPX5My8eWk3cOsG5rZ4nxId6YGO6Nbh7899vcsAjpCIsQEZHxuFZUga1JtwdYFyj/HGAd5C3DxHA5xgR5wYkDrM0Ci5CWOFiaiMh41avUOJxehNgT2dh/vqBhgLWNpQQP9fTCpHA5wn3bcoC1CWMR0hGeESIiMm43y2uw41QOYk9kI72wvGF7J1d7TAqX49HeXOfMFLEI6QiLEBGRafhjgPWWE9nYfToXlbUqAICFRMCwAHdMjvDBgC6uXOfMRLAI6QiLEBGR6Smvqcee07nYlJiNlOyShu3ebW0xKUyOieFyuDvZiBeQtMYipCMsQkREpu1CvhKbE7OxPfk6lL+vcyaVCHigWztMifTBID83niUyQixCOsIiRERkHqrrVNibdnuds8Rrf96G793WFpMjfPB4mDfaOfIskbFgEdIS7xojIjJf6QVl2JiYhe+S/jxLZCERMLyHO57q0wF9O7nwjjMDxyKkIzwjRERkvqrrVPjhdB42Hs9EclZJw/bObvZ4qk8HPNrbGzJbzktkiFiEdIRFiIiIAOB8nhIbjmdiR3IOKn6/48zWUoqxwV6Y2rcDenjJRE5If8UipCMsQkRE9Fdl1XXYeSoH649l4WJBWcP2sA5tMa2fL0b28ICVhUTEhASwCOkMixAREd2LRqPBiWu38N9jmfgxLa9h9mo3R2s8GeGDJyN9eAu+iFiEdIRFiIiI/k6hshobE7Ow4XgWbpTdXuPMQiJgdE9PTO/vi94+bUVOaH5YhLTEu8aIiKi5auvV+PlsPr5NuIYT1241bA+St8HT/X0xKtCTl81aCYuQjvCMEBERtcSZnFKs++0adqXkolalBgC0c7TGU306YEqkD1wcuL6ZPrEI6QiLEBERaaOovAabjmfhv8cyUfj7ZTMrCwkeDWmPpwd0hL+7o8gJTROLkI6wCBERkS7U1qvx45k8fHUkA6evlzZsH+jniqcHdMRgLuWhUyxCOsIiREREuqTRaHAy8xa+OpyBfefy8fvNZujSzgEzB3bE2OD2sLGUihvSBLAI6QiLEBER6Ut2cSXW/XYNsSeyUV5zeykPVwcrRPf1xVN9OqCtvZXICY0Xi5COsAgREZG+KavrEJuYja+PZiC3tBoAYGMpweOhcjwzsCM6uNiLnND4sAjpCIsQERG1ljqVGnvT8rA6/irO5ioBABIBGBXoiWcHdUKQvI24AY0Ii5CWOI8QERGJRaPRIOHqTayOv4qDF280bO/TyRnPDe6MIf5uEAQOrG4Mi5CO8IwQERGJ6UK+Eqvjr2JXSm7DMh7dPBzx3OBOeLiXFyylnKDxXliEdIRFiIiIDEFuSRXWHsnApsQsVNSqAADebW3x7KBOeDxUDlsr3mn2VyxCOsIiREREhqS0sg7rj2di7ZEM3KyoBQC42FthRn9fTO3jC5mdpcgJDQOLkI6wCBERkSGqrlNh68lsrIq/iuu3qgAA9lZSPNW3A54Z0Alujua9hAeLkI6wCBERkSGrV6nxw+k8rDx4BRcLygAA1hYSTAqX49lBneDd1k7khOJgEdIRFiEiIjIGarUGv14oxIoDl5GSXQIAsJAIGBvcHrOHdkZnNwdxA7YyFiEdYREiIiJj8set958fuIIjl4sAAIIAPNTTE3OGdkF3T/P4LGMR0hEWISIiMlYp2SVQHLiMX84VNGyL6u6OeQ90MfnJGVmEdIRFiIiIjN35PCUUBy5jT1oe/vjEH+Tvhhce6IIwX2dxw+kJi5COsAgREZGpuHKjHJ8fuIKdKTlQ/T45Y7/OLnjhQT/06eQicjrdYhHSEpfYICIiU5V1sxIrD13BtqRs1KluV4CIjs6Y/6Af+nV2MYnlO1iEdIRnhIiIyFTllFThi4NXEHsiG7UqNQCgt08bvBjlj4F+rkZdiFiEdIRFiIiITF1+aTW+OHQFmxKzUFN/uxCF/F6IBhlpIWIR0hEWISIiMheFymqsir+K9ccyGwpRsLwNXozyw2AjW/GeRUhHWISIiMjcFJZVY/Whq1h/PBPVdcZZiFiEdIRFiIiIzNWNshqsjr+C/x77sxAZyyUzFiEdYREiIiJzd6OsBqsOXbnjDFFvnzZYMMwfA7oYZiFiEdIRFiEiIqLbCsuqserQnWOIwn3bYsEwf/Tr7CpyujuxCOkIixAREdGdCpXVWHnoCjYcz0Lt74WoTydnxAzrioiOhjFTNYuQjrAIERER3Vt+aTU+P3gZmxP/nIdooJ8rFgzzR2+ftqJmYxHSERYhIiKixuWUVEFx4DK2nMhG/e9Ldwzt6oaYYV3R01smSiYWIR1hESIiImqa7OJKfPZrOr5L/nMtsxE93LFgmD+6ebTuZyiLkI6wCBERETVPRlEFlselY2dKDjQaQBCAh3t54cUoP3R2c2iVDCxCOsIiRERE1DLpBWVYtj8de9LyAAASAXi0tzfmP+gHubOdXt+bRUhHWISIiIi0cza3FJ/8cgn7zxcCACylAiaFyzHvAT+4O9no5T1ZhLSkUCigUCigUqlw6dIlFiEiIiItncq6haW/XMLh9CIAgLWFBFP7dMCsIZ3h4mCt0/diEdIRnhEiIiLSrWNXb+LjfRdx4totAED7NraIf3UopBLdzVDdnM9vC529KxEREdHf6NPJBVue64v49CJ8vO8iRvf01GkJai4WISIiImpVgiBgsL8bBvm5Nsw9JBYWISIiIhKFIAiwlIq7aKtE1HcnIiIiEhGLEBEREZktFiEiIiIyWyxCREREZLZYhIiIiMhssQgRERGR2WIRIiIiIrPFIkRERERmi0WIiIiIzBaLEBEREZktFiEiIiIyWyxCREREZLZYhIiIiMhscfX5Rmg0GgCAUqkUOQkRERE11R+f2398jjeGRagRZWVlAAC5XC5yEiIiImqusrIyyGSyRp8jaJpSl8yUWq1Gbm4uHB0dIQhCw/bw8HCcOHHinq+532P/u12pVEIulyM7OxtOTk66D98Mjf08rb2/5ry2Kc9tybFq7DEeR92/lsfxTro+jtrsszWPY2OPG+NxBAznd9Icj6NGo0FZWRm8vLwgkTQ+CohnhBohkUjg7e1913apVHrfg3O/x+633cnJSfRf2MZ+ntbeX3Ne25TntuRYNfYYj6PuX8vjeCddH0dt9tmax7Gxx43xOAKG8ztprsfx784E/YGDpVtgzpw5zX6ssdeITdfZtNlfc17blOe25Fg19hiPo+5fy+N4J31ka+k+W/M4Nva4MR5HwHB+J3kcG8dLYyJRKpWQyWQoLS01iL9cqGV4HE0Dj6Np4HE0Da19HHlGSCTW1tZYvHgxrK2txY5CWuBxNA08jqaBx9E0tPZx5BkhIiIiMls8I0RERERmi0WIiIiIzBaLEBEREZktFiEiIiIyWyxCBi47OxtDhgxBQEAAevXqha1bt4odiVpo/PjxaNu2LSZMmCB2FGqGH374AV27doWfnx/WrFkjdhzSAn8HjZ8+PhN515iBy8vLQ0FBAYKDg5Gfn4/Q0FBcunQJ9vb2YkejZjp48CDKysrwzTffYNu2bWLHoSaor69HQEAADhw4AJlMhtDQUPz2229wcXEROxq1AH8HjZ8+PhN5RsjAeXp6Ijg4GADg4eEBV1dXFBcXixuKWmTIkCFwdHQUOwY1Q2JiInr06IH27dvDwcEBo0aNwr59+8SORS3E30Hjp4/PRBYhLcXHx2PMmDHw8vKCIAjYuXPnXc9RKBTw9fWFjY0NIiMjkZiY2KL3SkpKgkqlglwu1zI1/a/WPI7UerQ9rrm5uWjfvn3D9+3bt0dOTk5rRKf/wd9R06DL46irz0QWIS1VVFQgKCgICoXino/HxsYiJiYGixcvRnJyMoKCgjBixAgUFhY2PCc4OBiBgYF3feXm5jY8p7i4GNOmTcPq1av1/jOZo9Y6jtS6dHFcyTDwWJoGXR1HnX4makhnAGh27Nhxx7aIiAjNnDlzGr5XqVQaLy8vzZIlS5q83+rqas3AgQM13377ra6iUiP0dRw1Go3mwIEDmscee0wXMamZWnJcjx49qhk3blzD4/Pnz9ds2LChVfLS/WnzO8rfQcPR0uOo689EnhHSo9raWiQlJSEqKqphm0QiQVRUFBISEpq0D41Gg+nTp+OBBx7A1KlT9RWVGqGL40iGpynHNSIiAmfOnEFOTg7Ky8vx448/YsSIEWJFpvvg76hpaMpx1MdnIouQHhUVFUGlUsHd3f2O7e7u7sjPz2/SPo4ePYrY2Fjs3LkTwcHBCA4ORlpamj7i0n3o4jgCQFRUFB5//HHs3bsX3t7e/AdaZE05rhYWFvj4448xdOhQBAcH46WXXuIdYwaoqb+j/B00bE05jvr4TLTQ6tWkdwMGDIBarRY7BunA/v37xY5ALfDII4/gkUceETsG6QB/B42fPj4TeUZIj1xdXSGVSlFQUHDH9oKCAnh4eIiUipqLx9E08biaDh5L0yDWcWQR0iMrKyuEhoYiLi6uYZtarUZcXBz69u0rYjJqDh5H08Tjajp4LE2DWMeRl8a0VF5ejsuXLzd8n5GRgZSUFDg7O8PHxwcxMTGIjo5GWFgYIiIisGzZMlRUVGDGjBkipqb/xeNomnhcTQePpWkwyOOok3vPzNiBAwc0AO76io6ObnjOZ599pvHx8dFYWVlpIiIiNMeOHRMvMN0Tj6Np4nE1HTyWpsEQjyPXGiMiIiKzxTFCREREZLZYhIiIiMhssQgRERGR2WIRIiIiIrPFIkRERERmi0WIiIiIzBaLEBEREZktFiEiIiIyWyxCREREZLZYhIiIiMhssQgREenA+PHj0bZtW0yYMEHsKETUDCxCREQ6MH/+fHz77bdixyCiZmIRIqJWM2TIELz44ovNer4gCBAEASkpKXrLpQtDhgyBo6PjPR97+eWXMW7cuPu+dvr06Q0/586dO/UTkIjuiUWIiAzazJkzkZeXh8DAwIZtgwcPhiAIePfdd+94rkajQWRkJARBwJtvvtnaUe8rJSUFvXr1uu/jn376KfLy8loxERH9wULsAEREjbGzs4OHh0fD9xqNBqdOnUKHDh2QlpZ2x3O/+eYb5ObmAgB69+6tswzBwcGor6+/a/u+ffvg5eX1t69PTU3FrFmz7vu4TCaDTCbTKiMRtQzPCBGRaPbs2QOZTIYNGzY0+TXp6ekoKytDdHT0HUWorKwMixYtwvTp0wEAoaGhAID8/HwIgoBPP/0UISEhsLGxQY8ePXDkyJE79puVlYXo6Gi4u7vD1tYWQUFBDc9JSUnBmTNn7vpqSgm6fv06ioqKAADDhg2DnZ0dunbtiuPHjzf5ZyYi/WERIiJRbNy4EZMnT8aGDRswZcqUJr8uKSkJdnZ2mDx5Mi5evIja2loAwFtvvYWwsDC4ubnBw8MDnp6eANAwtmjt2rVYtmwZUlJS4OPjgylTpkCtVgMAMjMzERERgaqqKuzatQunT5/G3Llz4eTkpPXP+cf7KxQKvPbaa0hNTYWPjw8WLlyo9b6JSHu8NEZErU6hUOCf//wndu/ejcGDBzfrtcnJyejVqxe6du0KGxsbXLhwAba2tli5ciWSk5Pxzjvv3HFZLDU1FZaWlvj+++/h6+sLAHj77bcRFhaGnJwcyOVyzJo1C3369MGWLVsaXufn59esXFFRUUhNTUVFRQW8vb2xdetW9O3bFykpKXB2dsaWLVvg6uoKAHjkkUewatWqZu2fiPSDRYiIWtW2bdtQWFiIo0ePIjw8vNmvT05ORu/evSEIAnr16oW0tDRs2rQJs2bNgp+fH5KSkjB+/PiG56ekpODRRx9tKEEA7jjTk5mZiR9//BGnTp3S6ufav3//PbenpKRg7NixDSUIADIyMtClSxet3o+IdIOXxoioVYWEhMDNzQ1r166FRqNp9uv/KELA7UHMy5Ytw8mTJ/H666+juroaFy5cuOOMUEpKCoKDg+/YR0JCAlxdXdG+fXukpKTAysrqrufoSkpKCvr06XPXNn29HxE1D4sQEbWqzp0748CBA/j+++8xb968Zr326tWrKCkpaSg6ISEhOHnyJJYsWQJHR0ekpqaivr6+YaB0VVUV0tPToVKpGvahVquxbNkyREdHQyKRwNLSEvX19aisrNTdD/m7srIyXL16FSEhIXdsZxEiMhwsQkTU6vz9/XHgwAF89913zZpgMSkpCVZWVg1zCkVHR+PGjRsNd4olJyfDzc0NcrkcAJCWlgZBELB+/XokJCTg/PnzmDRpEkpKSvCvf/0LABAZGQmZTIZZs2bh/PnzOHfuHL744gukp6dr/XOmpqZCKpWiZ8+eDdsyMzNx69YtFiEiA8EiRESi6Nq1K3799Vds2rQJL730UpNek5ycjMDAQFhaWgIALC0t4erqCkEQGh7/69mXlJQUdOvWDa+99hoee+wxhIWFQaVS4dChQ2jTpg0AwMXFBbt370Z6ejrCw8MxYMAA7Nq1C+3atdP6Z0xJSWkY1P2HU6dOoU2bNneMWSIi8QiallykJyJqBUOGDGkYB9QSc+bMwa1bt7Bx40bdBtMTQRCwY8eORpfjICLd4hkhIjJon3/+ORwcHO6aRbop/m5pC0Px/PPPw8HBQewYRGaJZ4SIyGDl5OSgqqoKAODj4wMrK6smv1aj0UAmk2Hz5s0YPXq0viLqRGFhIZRKJQDA09MT9vb2IiciMh8sQkRERGS2eGmMiIiIzBaLEBEREZktFiEiIiIyWyxCREREZLZYhIiIiMhssQgRERGR2WIRIiIiIrPFIkRERERmi0WIiIiIzBaLEBEREZktFiEiIiIyWyxCREREZLZYhIiIiMhs/T9Vf6nM6lheJwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(hm.k_hm, hm.power_auto_tracer)\n", "plt.xscale(\"log\")\n", "plt.yscale(\"log\")\n", "plt.xlabel(\"k [$Mpc^{-1} h$]\")\n", "plt.ylabel(r\"$\\rm P(k) \\ [{\\rm Mpc^3}h^{-3}]$\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And in temperature units:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAHACAYAAABUC+fAAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVExJREFUeJzt3XtcVHX+P/DXzADDRUARAbmJd8MLozAgW5YWRlZes3XNFKHMXLIL1a625eZu5u52WXZzii4iXtK8pVZu/izS8Bo3By8oSpICyggit+E6M+f3B9t81xAbYIYzM7yejwd/nDNnznlPJ5wXn8/nfD4SQRAEEBEREdGvkopdABEREZGtYHAiIiIiMhGDExEREZGJGJyIiIiITMTgRERERGQiBiciIiIiEzE4EREREZmIwYmIiIjIRAxORERERCZicCIiIiIyEYMTERERkYkcxC7A2hQXF2P+/Pm4du0aHBwc8Nprr+HRRx81+f0GgwFXrlyBu7s7JBKJBSslIiIicxEEAbW1tfD394dU2n67koSL/N7s6tWr0Gg0UCgUKCsrQ3h4OM6fPw83NzeT3l9SUoKgoCALV0lERESWUFxcjMDAwHZfZ4vTL/Tv3x/9+/cHAPj5+cHb2xuVlZUmByd3d3cArf/hPTw8LFYnERERmU9NTQ2CgoKM3+PtsbnglJGRgbfeegs5OTm4evUqdu3ahRkzZtx0jEqlwltvvYWysjKEhYXhvffeQ2RkZIevlZOTA71e36EWpJ+75zw8PBiciIiIbMyvDbOxueCk1WoRFhaGhIQEzJo1q83rW7duRVJSElJSUhAVFYXk5GTExsaioKAAPj4+AACFQgGdTtfmvfv374e/vz8AoLKyEgsWLMDHH39823qamprQ1NRk3K6pqenKxyMiIiIrZtNjnCQSSZsWp6ioKCiVSqxZswZA62DtoKAgLF26FMuWLTPpvE1NTZg8eTIWLVqE+fPn3/bY119/HStXrmyzv7q6mi1ORERENqKmpgaenp6/+v1tV9MRNDc3IycnBzExMcZ9UqkUMTExOHbsmEnnEAQBCxcuxL333vuroQkAli9fjurqauNPcXFxp+snIiIi62ZXwamiogJ6vR6+vr437ff19UVZWZlJ5zhy5Ai2bt2K3bt3Q6FQQKFQ4NSpU+0eL5fL4eHhgY0bN2L8+PG47777uvQZiIiIyHrZ3BgnS7vrrrtgMBg6/L7ExEQkJiYam/qIiIjI/thVi5O3tzdkMhk0Gs1N+zUaDfz8/Cx6bZVKhdDQUCiVSoteh4iIiMRjV8HJyckJ4eHhSE9PN+4zGAxIT09HdHS0Ra+dmJiI/Px8ZGVlWfQ6REREJB6b66qrq6tDYWGhcbuoqAhqtRpeXl4IDg5GUlIS4uLiEBERgcjISCQnJ0Or1SI+Pl7EqomIiMge2Fxwys7OxqRJk4zbSUlJAIC4uDikpaVhzpw5KC8vx4oVK1BWVgaFQoF9+/a1GTBubiqVCiqVCnq93qLXISIiIvHY9DxO1sjUeSCIiIjIevTIeZzExMHhRERE9o8tTmbGFiciIiLbwxYnIiIisivFlfX45NBFUWuwucHhRERE1POoi6vw5PosVNQ1w8PZEb9VBolSB1uczIRjnIiIiCzj61NXMefDY6ioa8Yd/T0wYZi3aLVwjJOZcYwTERGReQiCgI8yLmL11+cAAJOG98N7j41DL7n5O8xM/f5mVx0RERFZnRa9ASv2nMGWzMsAgAXRA7Di4VA4yMTtLGNwIiIiIqtS09iCxE9zcehCBSQS4LWHQhF/ZwgkEonYpTE4mQtnDiciIuq6khv1SEjLwnlNHVwcZfj33LGYHGrZ1T86gmOczIxjnIiIiDonr7gKT6zPRkVdE3zc5Vgbp8ToQM9uuTbHOBEREZHN+H9nyvDcZyfQ2GLACD93pC5Uwr+3i9hltcHgRERERKIRBAFrDxdh1X/OQhCAe4b1w5rHxsLd2VHs0m6JwYmIiIhEodMbsPLLfGw8fgkAMC8qGCunjRT9ybnbYXAyEw4OJyIiMl1dkw5LN+fiQEE5JBLglSl34MkJA63iybnb4eBwM+PgcCIiotsrq25EQloW8q/WwNlRiuQ5Cjwwqr+oNXFwOBEREVmd/Cs1SEjLQllNI7x7OeGTOCUUQb3FLstkDE5ERETULQ6cu4ZnNudC26zHEJ9eWLdQiSAvV7HL6hAGJyIiIrK4jccv4c97TsMgAL8Z3BcfPB4OTxfrfHLudhiciIiIyGIMBgGrvz6Ljw8VAQBmhwfizZmj4eRgvU/O3Q6DExEREVlEQ7MeL2xVY9+ZMgDAi5OH4Zl7h1j9k3O3w+BEREREZldR14Qn12dDXVwFJ5kUbz06BtMVAWKX1WUMTmbCeZyIiIhaFV6rQ3xaJoorG+Dp4oiP5ocjalBfscsyC87jZGacx4mIiHqy4xev46kN2ahp1CHYyxXr4pUY3K+X2GX9Ks7jRERERN1q94lSvLwjDy16AeOCe+PjBRHo20sudllmxeBEREREXSIIAtZ8V4h3vjkPAHhwtB/e/a0Czo4ykSszPwYnIiIi6rQWvQF/2nUK27JLAACL7x6EPz4wAlKp7T45dzsMTkRERNQpNY0t+P2mXBwurIBUAqycPgrzxw8QuyyLYnAiIiKiDrtS1YD4dVko0NTC1UkG1WPjMGmEj9hlWZxtTttpQVVVVYiIiIBCocCoUaPw8ccfi10SERGRVTldWo0ZqiMo0NTCx12ObYuje0RoAtji1Ia7uzsyMjLg6uoKrVaLUaNGYdasWejb1z7mnyAiIuqKAwXX8MynrQv1Dvd1R2q8EgG9XcQuq9swOP2CTCaDq2vrSs1NTU0QBAGc6oqIiAjY/MNlvLbnNPQGAXcOaV2o18PZ9hbq7Qqb66rLyMjA1KlT4e/vD4lEgt27d7c5RqVSISQkBM7OzoiKikJmZmaHrlFVVYWwsDAEBgbi5Zdfhre3t5mqJyIisj0Gg4C/7zuHV3adgt4g4JFxgVi3MLLHhSbABoOTVqtFWFgYVCrVLV/funUrkpKS8Oc//xm5ubkICwtDbGwsrl27Zjzm5/FLv/y5cuUKAKB3797Iy8tDUVERNm/eDI1G0y2fjYiIyNo06fR4fqsaHxz8EQDwfMxQvP3oGDg52FyEMAubXnJFIpFg165dmDFjhnFfVFQUlEol1qxZAwAwGAwICgrC0qVLsWzZsg5f4/e//z3uvfdezJ49+5avNzU1oampybhdU1ODoKAgLrlCREQ2r6q+GU9tzEFmUSUcpBL87ZExmB0eKHZZFmHqkit2FRebm5uRk5ODmJgY4z6pVIqYmBgcO3bMpHNoNBrU1tYCAKqrq5GRkYHhw4e3e/zq1avh6elp/AkKCurahyAiIrICxZX1eOSDo8gsqoS73AFp8ZF2G5o6wq6CU0VFBfR6PXx9fW/a7+vri7KyMpPOcenSJUyYMAFhYWGYMGECli5ditGjR7d7/PLly1FdXW38KS4u7tJnICIiEtvJkirMfP8ofizXor+nM7YvicZdQzneF+BTdW1ERkZCrVabfLxcLodcLodKpYJKpYJer7dccURERBaWflaDZzafQEOLHnf098C6hUr4eTqLXZbVsKsWJ29vb8hksjaDuTUaDfz8/Cx67cTEROTn5yMrK8ui1yEiIrKUjccvYdGGbDS06HH3sH7Ytng8Q9Mv2FVwcnJyQnh4ONLT0437DAYD0tPTER0dbdFrq1QqhIaGQqlUWvQ6RERE5mYwCPjb1+fw2u7TMAjAbyMCsTYuAu49cLqBX2NzXXV1dXUoLCw0bhcVFUGtVsPLywvBwcFISkpCXFwcIiIiEBkZieTkZGi1WsTHx1u0rsTERCQmJhpH5RMREdmCJp0eL28/iS/yWqfkeSFmGJ69bwgkEonIlVknmwtO2dnZmDRpknE7KSkJABAXF4e0tDTMmTMH5eXlWLFiBcrKyqBQKLBv3742A8bNjWOciIjI1lTXt+Cpjdn4oQdMN2AuNj2PkzUydR4IIiIiMZXcqEf8uixcuFaHXnIHfPD4OEwY2k/sskRj6ve3zbU4ERERUdecLq1GfFoWymub4Oshx7qFkQj15x/7pmBwMhN21RERkS3IOF+OJZtyoG3WY7ivO9bFK+Hf20XssmwGu+rMjF11RERkrbZnF2P556egMwiIHtQXKfPD4enCJ+cAdtURERHRfwmCgPe+K8S735wHAExX+OMfs8dA7iATuTLbw+BERERkx3R6A17bcxpbMluXBFsycTBevn84pFJON9AZDE5mwjFORERkbeqbdXhm8wl8d+4apBJg5bSRmB8dInZZNo1jnMyMY5yIiMgaVNQ1ISEtCydLquHsKMW/fzcW94+07PJjtoxjnIiIiHqoogot4lIzcbmyHl5uTvgkLgLjgvuIXZZdYHAiIiKyIycu38AT67NRqW1GsJcr1idEYqC3m9hl2Q27WuRXTFzkl4iIxPZtvgZzPz6OSm0zxgR6YueS3zA0mRnHOJkZxzgREZEYPv3hEl7bfRoGAZg4vB9Uj42Dm5wdS6biGCciIqIeQBAEvPvNebz3XSEA4LcRgVg1czQcZexUsgQGJyIiIhvVojfglc9PYXtOCQDg2fuG4oWYoZBIOEeTpTA4ERER2SBtkw6Jm3NxsKAcUgmwauZozI0MFrssu8fgZCacAJOIiLrLL+doWjN3HGJCfcUuq0fg4HAz4+BwIiKypEvXtViQmolL1+vRx9URaxcqOUeTGXBwOBERkZ05WVKF+HVZuK5tRmAfF6xPiMTgfr3ELqtHYXAiIiKyAd+fL8eSTTmob9YjtL8H0uKV8PFwFrusHofBiYiIyMp9nluCP+w4CZ1BwJ1D+iLl8XC4OzuKXVaPxOBERERkpQRBwEcZF7H663MAgGlh/nj70TA4OXCOJrEwOBEREVkhg0HAG3vPIvVIEQBg0YSBWD7lDkilnKNJTAxOREREVqZJp0fStjzsPXkVAPDqQ3fgyQmDRK6KAAYns+E8TkREZA61jS1YvDEHR3+8DkeZBG8/GobpigCxy6L/4jxOZsZ5nIiIqLOu1TYifl0WzlypgZuTDCnzwzFhaD+xy+oROI8TERGRDSmq0GJB6g8ormyAdy8nrFsYidGBnmKXRb/A4ERERCSy/53YMtjLFRufiMSAvm5il0W3wOBEREQkoozz5Xj6vxNbjgrwwLqFkejnLhe7LGoHgxMREZFI9qhL8eK2POPElh/Oj0AvOb+arRnvDhERkQhSDxfhL1/lAwAeHtMf7/w2DHIHmchV0a/h1KPtqK+vx4ABA/DSSy+JXQoREdkRQRDw933njKFp4W9C8O/fjWVoshFscWrHqlWrMH78eLHLICIiO6LTG7D881PYnlMCAHg5djh+P3EwJBLOBm4r2OJ0CxcuXMC5c+cwZcoUsUshIiI70dCsx+KNOdieUwKpBPjHI2OQOGkIQ5ONsbnglJGRgalTp8Lf3x8SiQS7d+9uc4xKpUJISAicnZ0RFRWFzMzMDl3jpZdewurVq81UMRER9XTV9S2Yv/YHpJ+7BrmDFB/Oj8BvlUFil0WdYHNddVqtFmFhYUhISMCsWbPavL5161YkJSUhJSUFUVFRSE5ORmxsLAoKCuDj4wMAUCgU0Ol0bd67f/9+ZGVlYdiwYRg2bBiOHj36q/U0NTWhqanJuF1TU9OFT0dERPamrLoRC1J/wHlNHTycHbB2oRLKEC+xy6JOsuklVyQSCXbt2oUZM2YY90VFRUGpVGLNmjUAAIPBgKCgICxduhTLli371XMuX74cmzZtgkwmQ11dHVpaWvDiiy9ixYoVtzz+9ddfx8qVK9vs55IrRET0Y3kdFqzNRGlVA3w95FifEIkRfvxusEamLrliV8GpubkZrq6u2LFjx01hKi4uDlVVVdizZ0+Hzp+WlobTp0/j7bffbveYW7U4BQUFMTgREfVwecVViE/LQqW2GYO83bA+IRJBXq5il0Xt6JFr1VVUVECv18PX1/em/b6+vjh37pxFrimXyyGXy6FSqaBSqaDX6y1yHSIish2HLpRj8cbW2cDHBHpi3UIl+vbibOD2wK6Ck7ktXLjQ5GMTExORmJhoTKxERNQzfXXyCl7YqkaLXsBdQ7yRMj+cs4HbEbu6k97e3pDJZNBoNDft12g08PPzE6kqIiLqKTYe+wkrvjgDQQAeGt0f787hbOD2xuamI7gdJycnhIeHIz093bjPYDAgPT0d0dHRFr22SqVCaGgolEqlRa9DRETWRxAEJH97Hq/taQ1N86KC8e+5nA3cHtlci1NdXR0KCwuN20VFRVCr1fDy8kJwcDCSkpIQFxeHiIgIREZGIjk5GVqtFvHx8Rati111REQ9k8EgYOWXZ7D+2CUAwHP3DcXzMUM5saWdsrnglJ2djUmTJhm3k5KSALQ+OZeWloY5c+agvLwcK1asQFlZGRQKBfbt29dmwLi5cXA4EVHP06wz4MXtefgy7wokEuD1qSMR95sQscsiC7Lp6QiskamPMxIRkW2rb9bh6U25yDhfDgepBO/8NgzTFQFil0Wd1COnIyAiIuoOVfXNiE/LwonLVXBxlOGDx8dh4nAfscuibsDgZCbsqiMi6hn+dwkVTxdHpC5UInxAH7HLom7CrjozY1cdEZH9ulheh/n/s4TKxieiMMzXXeyyyAzYVUdERGRGp0urEZeaievaZgz0dsMGLqHSI9nVPE5i4jxORET269iP1/G7j47jurYZI/09sP3paIamHopddWbGrjoiIvuy/0wZntlyAs06A8YP8sLHCyLg7uwodllkZuyqIyIi6qIdOSX4486T0BsETA71xXtzx8LZkbOB92QMTkRERLfwyaGLeGPvWQDA7PBA/G3WaDjIOMKlp2NwMhNOR0BEZB8EQcDb+wugOvAjAODJuwbilQfvgFTKJVSIY5zMjmOciIhsl94g4LU9p7H5h8sAgJdjh+P3Ewdz3bkegGOciIiIOqBZZ8AL29TYe/IqJBLgjRmjMC9qgNhlkZVhcCIioh6vvlmHxRtzcOhCBRxlEvxzjgIPj/EXuyyyQgxORETUo1XVNyMhLQu5/113LmV+OO4Z1k/ssshKMTgREVGPpalpxIK1mSjQ1HLdOTIJn6s0E84cTkRkWy5d12J2ylEUaGrh4y7HtsXRDE30q/hUnZnxqToiIut39moNFqRmory2CcFervj0ySguodLD8ak6IiKiW8i5VIn4dVmoadRhhJ87NjwRCR93Z7HLIhvB4ERERD3GwYJreHpTDhpbDIgY0Adr45TwdOW6c2Q6BiciIuoRvsi7gqStaugMAiYO74cP5oXDxYnrzlHHMDgREZHd23T8El7bcxqCAEwN88c7j4bByYHPR1HHMTgREZHdEgQB7x/8EW/9vwIAwOPjg7Fy2ijIuO4cdRKDk5lwkV8iIusiCAJW7T2LTw4XAQCemTQEL94/jOvOUZdwOgIz43QERETi0+kNWPb5KezIKQEAvPZwKJ64a6DIVZE143QERETUIzW26PHslhPYn6+BTCrB3x8Zg9nhgWKXRXaCwYmIiOxGXZMOi9Zn49jF63BykEL12DhMDvUVuyyyIwxORERkFyq1zVi4LhMnS6rRS+6AjxdEIHpwX7HLIjvD4ERERDavtKoB89f+gIvlWni5OSEtXokxgb3FLovsEIMTERHZtMJrdZi/9gdcrW6Ev6czNj4ZhcH9eoldFtkpBiciIrJZJ0uqEJeaiRv1LRjczw0bn4iCf28XscsiO8bgdAshISHw8PCAVCpFnz59cODAAbFLIiKiXzhaWIFFG7KhbdZjTKAn0uIj4eXmJHZZZOcYnNpx9OhR9OrFpl4iImu07/RVPLtFjWa9Ab8Z3BcfLYhALzm/0sjy+H8ZERHZlM8yL+OVXadgEIDYkb741+/GwtmRi/VS97C5FQ4zMjIwdepU+Pv7QyKRYPfu3W2OUalUCAkJgbOzM6KiopCZmdmha0gkEtxzzz1QKpX49NNPzVQ5ERF1Reu6c4VY9nlraPqdMgjvzwtnaKJuZXMtTlqtFmFhYUhISMCsWbPavL5161YkJSUhJSUFUVFRSE5ORmxsLAoKCuDj4wMAUCgU0Ol0bd67f/9++Pv74/DhwwgICMDVq1cRExOD0aNHY8yYMbesp6mpCU1NTcbtmpoaM31SIiL6mcEg4M3//N+6c0smDsYfYodz3Tnqdja9Vp1EIsGuXbswY8YM476oqCgolUqsWbMGAGAwGBAUFISlS5di2bJlHb7Gyy+/jJEjR2LhwoW3fP3111/HypUr2+znWnVERObRojdg2c5T2Jnbuu7cqw/dgScnDBK5KrI3pq5VZ3NddbfT3NyMnJwcxMTEGPdJpVLExMTg2LFjJp1Dq9WitrYWAFBXV4fvvvsOI0eObPf45cuXo7q62vhTXFzctQ9BRERGjS16LNmUg525JZBJJXjn0TCGJhKVzXXV3U5FRQX0ej18fW9el8jX1xfnzp0z6RwajQYzZ84EAOj1eixatAhKpbLd4+VyOeRyOVQqFVQqFfR6fec/ABERGVXXt+DJDVnI+ukG5P9ddy6G686RyOwqOJnDoEGDkJeX1+H3JSYmIjEx0djUR0REnVdW3Yi41EwUaGrh7uyAtXFKRA70ErssIvsKTt7e3pDJZNBoNDft12g08PPzs+i12eJERGQeP5bXYcHaTJRWNcDHXY4NT0RihB/HjJJ1sKsxTk5OTggPD0d6erpxn8FgQHp6OqKjoy167cTEROTn5yMrK8ui1yEismd5xVV4NOUYSqsaMNDbDTuX/IahiayKzbU41dXVobCw0LhdVFQEtVoNLy8vBAcHIykpCXFxcYiIiEBkZCSSk5Oh1WoRHx8vYtVERPRrDl0ox+KNOahv1mN0gCfWxSvh3UsudllEN7G54JSdnY1JkyYZt5OSkgAAcXFxSEtLw5w5c1BeXo4VK1agrKwMCoUC+/btazNg3NzYVUdE1Hm7T5Tipe150BkE3DXEGynzw7mEClklm57HyRqZOg8EERG1+jjjIlb95ywA4OEx/fHOb8Mgd+Bs4NS9TP3+ZpwnIiJR/HI28IQ7B+LVh+6AVMrZwMl6MTiZCbvqiIhM16wz4OUdedijvgIAWD5lBJ66exCXUCGrx646M2NXHRHR7dU16bBkUw4OXaiAg1SCf8weg1njAsUui3o4dtUREZHVuVbTiIXrspB/tQauTjJ88Hg47hnWT+yyiEzG4GQm7KojIrq9wmu1iEvNQmlVA/q6OSF1oRJhQb3FLouoQ9hVZ2bsqiMiaivrp0o8uT4b1Q0tGOjthrR4JQb0dRO7LCIjdtUREZFV+PrUVTy3VY1mnQFjg3tjbZwSXm5OYpdF1CkMTkREZDGph4vw1735EARgcqgv/v27sXBx4hxNZLsYnMyEY5yIiP6P3iBg1d6zSD3SOkfT/PED8Pq0kZBxjiaycRzjZGYc40REPV19sw7PfabGN/kaAMAfHxiBp+/hHE1k3Sw6xumLL77o8HsmT54MFxeXzlyOiIhsxLXaRjy5PhsnS6rh5CDFO4+GYWqYv9hlEZlNp4LTjBkzOnS8RCLBhQsXMGjQoM5cjoiIbMAFTS0WrmudbqCPqyM+XhCBiBAvscsiMqtOj3EqKyuDj4+PSce6u7t39jJERGQDjhZWYPGmHNQ26jDQ2w3rFioR4s3pBsj+dCo4xcXFdajb7fHHH+d4HyIiO/VZ5mW8uvs0dAYBEQP64OMFEejD6QbITnFwuJn871N158+f5+BwIrJ7eoOAv319Fh8fan1ybmqYP96aPQbOjpxugGyPqYPDGZzMjE/VEVFPoG1qfXLu27OtT869EDMMz943hE/Okc0y9ftb2tETNzQ0oLS0tM3+M2fOdPRURERkg65UNWB2yjF8e1YDJwcp/j13LJ6LGcrQRD1Ch4LTjh07MHToUDz00EMYM2YMfvjhB+Nr8+fPN3txRERkXdTFVZiuOoKzV2vg3UuOz54aj2mcboB6kA4FpzfeeAM5OTlQq9VYt24dnnjiCWzevBkAwB4/IiL79nluCX774TGU1zZhhJ879jxzJ8YF9xG7LKJu1aGn6lpaWuDr6wsACA8PR0ZGBmbOnInCwkI20RIR2Smd3oC/7ztnHAQ+OdQX/5yjQC85V+2inqdDLU4+Pj44efKkcdvLywvffPMNzp49e9N+IiKyD9X1LYhPyzKGpmfvHYIPHw9naKIeq0PBaePGjW0mvXRycsKWLVvw/ffftzm+rq6ua9XZEJVKhdDQUCiVSrFLISIyi8JrtZiuOoxDFyrg4iiD6rFxSLp/OKRcqJd6sE5PR/DPf/4TL7zwQruv19bW4oEHHsCRI0c6XZwt4nQERGQP9p8pQ9K2PNQ16RDQ2wUfLQjHSH9PscsishiLLvILAK+88gr69u2LBQsWtHlNq9XigQcewPXr1zt7eiIiEoHeIOCd/QV4/+CPAIDIgV74YN449O0lF7kyIuvQ6eC0ceNGzJ8/H71798a0adOM+7VaLWJjY1FeXn7L7jsiIrJOldpmPPfZCRy6UAEAiL8zBK88eAccZR2e8o/IbnU6OM2ePRtVVVWYO3cu9u7di4kTJxpbmjQaDb7//nv079/fnLUSEZGFnCqpxtObclBa1QAXRxn+9shoTFcEiF0WkdXp0mMRTz75JCorKzF9+nTs2bMHK1aswJUrV/D999/D358TohER2YJtWcV4dc9pNOsMCOnripT54RjhxzGaRLfS5edJ//CHP6CyshL33XcfQkJCcPDgQQQGBpqjNiIisqD6Zh1e3X0an+e2LqMVc4cP3vmtAp4ujiJXRmS9Oh2cZs2addO2o6MjvL298dxzz920//PPP+/sJYiIyEIKymqRuDkXhdfqIJUAL94/HEvuGcypBoh+RaeDk6fnzY+lzp07t8vFWIuioiIkJCRAo9FAJpPh+PHjcHNzE7ssIiKz2JZdjBV7TqOxxQBfDzn+/buxiBrUV+yyiGxCp+dxsmf33HMP3njjDUyYMAGVlZXw8PCAg4NpGZPzOBGRtapv1uG13WewM7cEADBhqDf+OUcBb041QGTy93ennjE9efIkDAaDycefOXMGOp2uM5fqdmfOnIGjoyMmTJgAoHVZGVNDExGRtTpVUo2H/30YO3NLIJUAL8cOx/r4SIYmog7qVHAaO3Zshya3jI6OxuXLlztzqTYyMjIwdepU+Pv7QyKRYPfu3W2OUalUCAkJgbOzM6KiopCZmWny+S9cuIBevXph6tSpGDduHN58802z1E1EJAa9QcAHB3/EzPeP4GKFFn4eztiyaDwSJw3heCaiTuhUU4ogCHjttdfg6upq0vHNzc2ducwtabVahIWFISEhoc0AdQDYunUrkpKSkJKSgqioKCQnJyM2NhYFBQXGdfYUCsUtW8D2798PnU6HQ4cOQa1Ww8fHBw888ACUSiUmT558y3qamprQ1NRk3K6pqTHTJyUi6porVQ1I2qbG8YuVAIApo/zw5szR6OPmJHJlRLarU8Hp7rvvRkFBgcnHR0dHw8XFpTOXamPKlCmYMmVKu6+/++67WLRoEeLj4wEAKSkp2Lt3L1JTU7Fs2TIAgFqtbvf9AQEBiIiIQFBQEADgwQcfhFqtbjc4rV69GitXruzkpyEisoz/nLqK5Z+fQnVDC1ydZHh96kg8GhEIiYStTERd0angdPDgQTOXYR7Nzc3IycnB8uXLjfukUiliYmJw7Ngxk86hVCpx7do13LhxA56ensjIyMDixYvbPX758uVISkoybtfU1BhDFxFRd6uqb8bKL/Ox60Tr3ExhgZ5I/t1YDPTmk8FE5mBXo54rKiqg1+vh6+t7035fX1+cO3fOpHM4ODjgzTffxN133w1BEHD//ffj4Ycfbvd4uVwOuVwOlUoFlUoFvV7fpc9ARNRZ6Wc1WP75KVyrbYJUAiyZOBjPxwzjWnNEZmRXwclcfq078FYSExORmJhofJyRiKi7VDe04K9f5WNHTus0A4P6ueHtR8MwLriPyJUR2R+7Ck7e3t6QyWTQaDQ37ddoNPDz87PotdniRERiOFhwDct2nkJZTSMkEmDRhEFImjwMzo4ysUsjskt21X7r5OSE8PBwpKenG/cZDAakp6cjOjraotdOTExEfn4+srKyLHodIiIAqKhrwnOfncDCdVkoq2nEQG83bF8cjVcevIOhiciCbK7Fqa6uDoWFhcbtoqIiqNVqeHl5ITg4GElJSYiLi0NERAQiIyORnJwMrVZrfMqOiMiWCYKA7dklWPWfs6huaIFEAsT/ZiBejh0OFycGJiJLs7nglJ2djUmTJhm3f36iLS4uDmlpaZgzZw7Ky8uxYsUKlJWVQaFQYN++fW0GjJsbu+qIyNJ+LK/DK5+fwg9FrfMyhfb3wOpZoxEW1Fvcwoh6EK5VZ2Zcq46IzK2xRY8PDv6IDw7+iGa9AS6OMrwweSgS7hwIBz4xR2QWpn5/21yLk7ViixMRmZsgCPh/ZzT461f5KK1qAADcM6wf3pgxCkFepq3cQETmxRYnM2OLExGZQ+G1Oqz88gwOXagAAPT3dMafHroDD43uz9m/iSyALU5ERDaotrEF731XiNTDRdAZBDjJpHjq7kH4/aTBcHXiP9lEYjPLb2FLSwvKyspQX1+Pfv36wcvLyxynJSLqMVr0BmzJvIzkby+gUtu6MHrMHT547eFQDOjL5VKIrEWng1NtbS02bdqEzz77DJmZmWhuboYgCJBIJAgMDMT999+Pp556Ckql0pz1Wi2OcSKizvh5HNM/9p3DxQotgNaZv197OBSThvuIXB0R/VKnxji9++67WLVqFQYPHoypU6ciMjIS/v7+cHFxQWVlJU6fPo1Dhw5h9+7diIqKwnvvvYehQ4daon6rwzFORGSqE5dv4M3/nEXWTzcAAN69nPBczDD8ThnE9eWIupmp39+dCk5z587Fq6++ipEjR972uKamJqxbtw5OTk5ISEjo6GVsEoMTEf2ac2U1eGf/eXyT37o8lLOjFIsmDMLiewajl5zjmIjEYNHg9L+uXbsGHx82J/+MwYmI2nOxvA7//PYCvjp5BYIASCXAI+MC8eL9w+Hn6Sx2eUQ9Wrc9VTd79mwcOHAAMlnbqf51Oh0cHHrGX08c40RE7SmurMe/0y9gZ24JDP/9U/WhMf3xQswwDPHpJW5xRNQhXW5xmjZtGoKCgqBSqW7af/36dTzyyCM4ePBgV05vc9jiREQ/+6lCi/cPFuLz3FLo/puYYu7wQdLk4Qj1578PRNak21qcNmzYgMjISKSmphrHMZ09exYPP/wwRowY0dXTExHZnAuaWqgOFOKLvCvGFqYJQ72RNHkYxgb3Ebc4IuqSLgen3r17Y+fOnZg4cSJGjRqFGzduYM6cOXjiiSfw1ltvmaNGIiKbcLq0Gu8fLMTXp8vwc1v+vSN88My9QzCOgYnILnQqOM2aNQsKhcL4M3r0aKxZswYPPvggGhsb8d577yE+Pt7ctRIRWR1BEHC4sAIffn8RhwsrjPtjR/pi6b1DMSrAU8TqiMjcOhWcBg8ejEOHDmHNmjWoqKhAnz59EBYWBkEQ8Nhjj2HcuHFoaWmBo6Ojueu1WhwcTtSz6PQG7D11FR9+fxH5V2sAADKpBA+P6Y/fTxyC4X7uIldIRJbQ5cHhpaWlUKvVN/1cvHgRDg4OGDFiBPLy8sxVq03g4HAi+1bT2IJtWcVIO/oTSm40AABcHGWYowzCE3cNRJCXq8gVElFndNvg8ICAAAQEBOChhx4y7qurq4Nare5xoYmI7Nel61qsO/ITtmcXQ9vc2rLs5eaEhb8JwfzxA9DHzUnkComoO3Sqxeny5csIDg42+fjS0lIEBAR09DI2iS1ORPZDEAQcu3gd6478hG/PaowDvof49ELCnQMxc2wAXJzazmFHRLbH1O/vTi2GpFQqsXjxYmRlZbV7THV1NT7++GOMGjUKO3fu7MxliIhEUdekw8bjlxCbnIHHPv4B3+S3hqaJw/thQ0IkvnnhbjwWFczQRNQDdaqrLj8/H6tWrcLkyZPh7OyM8PBw+Pv7w9nZGTdu3EB+fj7OnDmDcePG4R//+AcefPBBc9dNRGR2hdfqsOn4JezIKUFdkw4A4Ookw8yxAYi/cyBn+Sairg0Ob2howN69e3H48GFcunQJDQ0N8Pb2xtixYxEbG4tRo0aZs1abwK46ItvSrDNgf34ZNv9wGUd/vG7cP8jbDfOjB+CR8EB4OPecJ4SJeqpuW+SXbsbgRGQbLl+vx+bMy9iRU4yKumYArYvu3neHLxZED8Cdg70hlUpErpKIuku3PVVHrTiPE5H1a9YZ8O1ZDbZkXsahC/83WaWPuxxzlEGYowxCYB9OJ0BE7etwi1NDQwMqKyvbPCV35swZjBw50qzF2SK2OBFZn/OaWmzNKsauE6Wo1La2LkkkwISh/fBYZDDuu8MHjrJOPStDRHbCIi1OO3bswPPPPw9vb28YDAZ8/PHHiIqKAgDMnz8fubm5XauaiMhMahtbsPfkVWzNLsaJy1XG/T7ucswOD8TvlMEI7svWJSLqmA4FpzfeeAM5OTnw9fVFTk4O4uLi8Morr+Cxxx4Dh0oRkdgMBgHHL17H9pwSfH36KhpbDAAAB6kE947wwRxlEO4Z1g8ObF0iok7qUHBqaWmBr68vACA8PBwZGRmYOXMmCgsLIZFwECURiePSdS125pRgZ24pSqsajPsH93PDoxFBmDUuAD7uziJWSET2okPBycfHBydPnsSYMWMAAF5eXvjmm28QFxeHkydPWqRAIqJbqW5owX9OXcXOnBJkX7ph3O/u7IBpYf6YHR4IRVBv/lFHRGbVocHhJSUlcHBwgJ+fX5vXjhw5gjvvvPOmfadPn+5xczlxcDiR5ej0BmRcKMfO3FJ8k69Bs661K04qAe4c4o1HI4Jwf6gvnB05ozcRdYxFBocHBga2+9rPoam2thZbtmzB2rVrkZOTA51O15FLiK6goABz5sy5aXvLli2YMWOGeEUR9WCCIOBUaTU+zy3FVyevGOdcAoBhvr3wyLhATFcEwM+TXXFEZHlmm8cpIyMDa9euxc6dO+Hq6ooJEyYgOzvbXKfvNsOHD4darQYA1NXVISQkBJMnTxa3KKIeqORGPfaor+Dz3BL8WK417u/r5oRpCn88Mi4QI/092BVHRN2qS8GprKwMaWlpWLt2La5evYrp06dj27ZtuP/++3Hu3Dns3r3bTGWK44svvsB9990HNzc3sUsh6hHqmnT4+tRV7MwtwfGLlcb9cgcpJof6Yta4AEwY2o9zLhGRaDodnKZOnYr09HRMmjQJr7/+OmbMmHFTwLDUX4EZGRl46623kJOTg6tXr2LXrl1tutFUKhXeeustlJWVISwsDO+99x4iIyM7fK1t27ZhwYIFZqqciG5FbxBw7Mfr2Jlbgn2ny9DQ0jr7vkQCRA30wqyxgXhgtB/XiyMiq9Dp4LR371489thjeP755xEREWHOmm5Lq9UiLCwMCQkJmDVrVpvXt27diqSkJKSkpCAqKgrJycmIjY1FQUEBfHx8AAAKheKWY6/2798Pf39/AK2DxI4ePYrPPvvMsh+IqIcqrqzH9pwS7MguxpXqRuP+gd5ueGRcAGaOC0RAbxcRKyQiaqvTweno0aNYu3Yt7r33XvTv3x/z5s3DvHnzMHjwYHPW18aUKVMwZcqUdl9/9913sWjRIsTHxwMAUlJSsHfvXqSmpmLZsmUAYBzDdDt79uzB/fffD2fn2w84bWpqQlNTk3G7pqbGhE9B1DM16fTYf0aDbdnFOFxYgZ+f6fVwdsDUMH88Eh6IsZxCgIisWKeD0/jx4zF+/HgkJydj69atSE1NxcqVK6FUKjFv3jxR1q1rbm5GTk4Oli9fbtwnlUoRExODY8eOdehc27Ztw1NPPfWrx61evRorV67scK1EPYUgCDhRXIU9J0rxRd4V3KhvMb5255C+mKMM5hQCRGQzuvxUnZubGxISEpCQkICCggKsXbsWb775JjQaTbf/1VhRUQG9Xm+c3fxnvr6+OHfunMnnqa6uRmZmJnbu3Pmrxy5fvhxJSUnG7ZqaGgQFBZleNJGd+rG8DntOlGJP3hVcul5v3O/n4YxHIwLxaHgQ14ojIptjtukIgNZH+f/xj39g9erV+PLLL5GammrO03cbT09PaDQak46Vy+WQy+VQqVRQqVTQ6/UWro7IejW26LFHXYpPf7iMkyXVxv0ujjLEjvTF9LEBuHtoP8ik7IojIttk1uD0M5lMhhkzZnT7pJHe3t6QyWRtQo9Go7nlbOfmlJiYiMTEROPMo0Q9yZWqBmw6fglbMi8bu+JkUgkmDPXGzLEBmBzqC1cni/xzQ0TUrezqXzInJyeEh4cjPT3dGNoMBgPS09PxzDPPWPTabHGinkYQBGT9dAPrj/6EfWfKoDe0jvQO6O2CBdED8Eh4ILx7yUWukojIvGwuONXV1aGwsNC4XVRUBLVaDS8vLwQHByMpKQlxcXGIiIhAZGQkkpOTodVqjU/ZWQpbnKinKK6sx64Tpfg8twQ//c/YpfGDvLDwNwMxOdSXXXFEZLdsLjhlZ2dj0qRJxu2fB2bHxcUhLS0Nc+bMQXl5OVasWIGysjIoFArs27evzYBxc2OLE9kzbZMO/7nFjN6uTjI8PKY/Fv5mIEL9uag1Edk/iSD8PJMKmYOpqysT2YK6Jh3WH/0JHx+6iKr/jl2SSIDoQX3xyLhAPDDKD25ym/v7i4ioDVO/v/kvHhG1oW3SYf2xn/BxxkXjYO8BfV3xaHggZ/Qmoh6NwclM2FVH9uCGthlbs4vxUcZFVGqbAbQugfLsfUMwLSyAY5eIqMdjV52ZsauObM3l6/XYn1+G/fkaZP9Uif8+HIeQvq5Yeu9QTFf4w0EmFbdIIiILY1cdEbWruqEFqYeLsO90GQo0tTe9dkd/DyTcGYKZYwMYmIiIfoHBiagHEQQBX+RdwV+/OouKutbFqWVSCaIGemFyqC9i7vBFkBeXQSEiag+Dk5lwjBNZu58qtHhtz2kculABABjUzw3PTBqC+0b4wtPVUeTqiIhsA8c4mRnHOJG1adLpkXLwIlQHC9GsM8DJQYqlk4bgqXsGQe4gE7s8IiKrwDFORD1cY4seO3JK8FHGRVyubJ3he8JQb/x1+iiEeLuJXB0RkW1icCKyMzWNLdh0/BJSD/9kHMfk3UuOFVNDMXVMf0gknFKAiKizGJzMhGOcSGyXr9fj08xL+PT4ZdQ16QC0Lri7aMJA/FYZBFcn/roTEXUVxziZGcc4UXcRBAFnrtRgf74G+8+U4VzZ/00rMMy3F5ZMHIyHx/jDkVMKEBH9Ko5xIrJT12oa8WHGRfy/M2UoudFg3C+TShA9qC/i7wzBpOE+kHKWbyIis2NwIrIRBoOAzZmX8fd951Db2NoV5+woxd1D+yF2pB/uu8MHvV2dRK6SiMi+MTgR2YDzmlos//wUci7dAACEBXpiycQhuGdYP7g4cUoBIqLuwuBkJhwcTpbQ2KKH6kAhUr7/ES16AW5OMrwUOxwLokO44C4RkQg4ONzMODiczOXMlWos3XICF8u1AICYO3zwl+mj4N/bReTKiIjsDweHE9koQRCw6fgl/HXvWTTrDPBxl2PltJF4YJQf52AiIhIZgxORFalpbMGynSfxn1NlAFpbmd6aHYY+bhz0TURkDRiciKzEyZIqJG7ORXFlAxykEiybMgJP3DWQrUxERFaEwYlIZA3NeqQeKULyt+fRohcQ2McFax4bB0VQb7FLIyKiX2BwIhJJs86ArdnFeC/9Aq7Vtq4p98BIP/x99hh4ujiKXB0REd0KgxNRN9MbBHyRV4p/fnMBlyvrAbSuKffi/cMwc2wAu+aIiKwYg5OZcB4n+jWNLXp8dfIqPsr4Eec1dQAA715yPHvfEMxRBkHuwIksiYisHedxMjPO40S/dLG8Dp/+cBk7ckpQ3dACAPBwdsDTEwdj4W9C4OrEv1+IiMTGeZyIRGQwCNifX4YNxy7h6I/XjfsDervgsahgPB41AJ6uHMdERGRrGJyIzEgQBHx/vhz/2FeA/Ks1AACJBLh3uA/mjQ/GPcN8uFQKEZENY3AiMpPcyzfw96/P4YeiSgCAu9wB86MH4LGoYAT2cRW5OiIiMgcGJ6IuKrxWh3/sO4f9+RoAgJODFAvGD8DvJw2BF2f8JiKyKwxORJ0kCAI2HLuEVf9pXVNOKgEeGReI5ycPQwAX4iUisksMTrfwz3/+E5988gkEQUBMTAz+9a9/cW4dukmlthl/2JGHb89eAwDcPawfXnvoDgz1dRe5MiIisiQGp18oLy/HmjVrcObMGTg6OuLuu+/G8ePHER0dLXZpZCWOFFbgha1qXKttgpODFK9MGYG434QwXBMR9QAMTreg0+nQ2NgIAGhpaYGPj4/IFZE1aNEb8O4355Hy/Y8QBGCITy+8N3cs7ujP+bqIiHoKqdgFdFRGRgamTp0Kf39/SCQS7N69u80xKpUKISEhcHZ2RlRUFDIzM00+f79+/fDSSy8hODgY/v7+iImJweDBg834CcgWVWqbMe+TH/DBwdbQNDcyGF8+cxdDExFRD2NzLU5arRZhYWFISEjArFmz2ry+detWJCUlISUlBVFRUUhOTkZsbCwKCgqMLUcKhQI6na7Ne/fv3w8XFxd89dVX+Omnn+Di4oIpU6YgIyMDd999t8U/G1mnC5paPLE+G5cr6+Eud8Bbj47BA6P6i10WERGJwOaC05QpUzBlypR2X3/33XexaNEixMfHAwBSUlKwd+9epKamYtmyZQAAtVrd7vu3b9+OIUOGwMvLCwDw0EMP4fjx4+0Gp6amJjQ1NRm3a2pqOvqRyIodLLiGpZtPoLZJh2AvV6QujMAQHw4AJyLqqWyuq+52mpubkZOTg5iYGOM+qVSKmJgYHDt2zKRzBAUF4ejRo2hsbIRer8fBgwcxfPjwdo9fvXo1PD09jT9BQUFd/hwkPkEQkHakCAlpWaht0iEyxAu7E+9kaCIi6uHsKjhVVFRAr9fD19f3pv2+vr4oKysz6Rzjx4/Hgw8+iLFjx2LMmDEYPHgwpk2b1u7xy5cvR3V1Nd5++20MHz4cQ4YM6dJnIPHp9Aa8tuc0Xv8yHwYBmB0eiI1PRnIySyIisr2uuu6watUqrFq1yqRj5XI55HI5XnzxRbz44ovG1ZXJNtU2tuD3n+bi0IUKSCTAsgdG4Km7B3GqASIiAmBnwcnb2xsymQwajeam/RqNBn5+fiJVRbbiSlUDEtKycK6sFi6OMvzrdwrcP5L/3xAR0f+xq646JycnhIeHIz093bjPYDAgPT3d4hNYqlQqhIaGQqlUWvQ6ZBmnSqoxQ3UE58pq0c9djm2LoxmaiIioDZtrcaqrq0NhYaFxu6ioCGq1Gl5eXggODkZSUhLi4uIQERGByMhIJCcnQ6vVGp+ys5TExEQkJiayq84GfZuvwdItJ9DQosdwX3ekxiu51hwREd2SzQWn7OxsTJo0ybidlJQEAIiLi0NaWhrmzJmD8vJyrFixAmVlZVAoFNi3b1+bAePmplKpoFKpoNfrLXodMq+0I0X4y1etg8AnDPWGat44eDg7il0WERFZKYkgCILYRdiTn1ucqqur4eHBWaWtVYvegJVfnsGm45cBAHMjg/CX6aPgKLOr3msiIjKRqd/fNtfiRNRVVfXNSNyciyOF1yGRAH98YAQW88k5IiIyAYOTmbCrzjb8WF6HJ9dno6hCC1cnGf71u7GYHGrZblwiIrIf7KozM3bVWa9DF8qR+Gkuahp1COjtgk/iIrhILxERAWBXHZGRIAjYcOwS/vJVPvQGAeED+iDl8XD0c5eLXRoREdkYBiczYVeddWps0ePV3aexI6cEADBrbADenDUazo4ykSsjIiJbxK46M2NXnfW4UtWApzfl4GRJNaQSYNmUEVg0gYPAiYioLXbVUY92/OJ1JH6ai+vaZvRxdcR7c8fhrqHeYpdFREQ2jsGJ7IogCFh/9Ce8sfcsdAYBof098OH8cAR5uYpdGhER2QEGJ7IbDc16/GnXKXx+ohQAMF3hj7/NGgMXJ45nIiIi82BwMhMODhdXUYUWSzbl4FxZLWRSCZZPGYEn7hrI8UxERGRWHBxuZhwc3v32nynDi9vyUNukg3cvOdY8NhbjB/UVuywiIrIhHBxOdk+nN+Dt/eeR8v2PAICIAX2gmjcOvh7OIldGRET2isGJbFJFXROe3XICR3+8DgBIuHMglj84gov0EhGRRTE4kc3JvXwDv9+Ui7KaRrg6yfCP2WPw8Bh/scsiIqIegMHJTDg43PJ+Xjrljb35aNELGNzPDSmPh2Oor7vYpRERUQ/BweFmxsHhllHfrMMrn5/CbvUVAMBDo/vj77PHoJec2Z+IiLqOg8PJbhRVaPH0xhwUaDjVABERiYvBiazad+c0eO4zNWobdejnLofqsXGIHOgldllERNRDMTiRVTIYBKgOFOLdb89DEFqnGnh/3jj4cKoBIiISEYMTWZ26Jh1e2paHfWfKAADzxw/Aaw+HwsmBUw0QEZG4GJzIqhRVaPHUhmxcuFYHJ5kUf50xEnOUwWKXRUREBIDByWw4HUHXfX++HM9szkVtow6+HnJ88Hg4xgX3EbssIiIiI05HYGacjqBz1h/9CSu/PAODAIQP6IMPOJ6JiIi6EacjIJug0xvwl6/yseHYJQDA7PBArJo5CnIHmciVERERtcXgRKKpaWzBM5tPION8OQDgjw+MwNP3DOL8TEREZLUYnEgUxZX1SEjLwoVrdXB2lCJ5jgIPjOovdllERES3xeBE3U5dXIUn0rJwXdsMXw85PlmgxOhAT7HLIiIi+lUMTtStDhRcw+835aKhRY+R/h5YG6eEnycHgRMRkW1gcKJusyOnBH/ceRJ6g4C7h/XDB/PGwY2L9BIRkQ3hVMy38Pbbb2PkyJEYNWoUNm3aJHY5Nk8QBHxw8Ee8tD0PeoOAmWMD8MmCCIYmIiKyOfzm+oVTp05h8+bNyMnJgSAImDRpEh5++GH07t1b7NJsksEg4K9787HuyE8AgMV3D8IfHxgBqZRPzhERke1hi9MvnD17FtHR0XB2doaLiwvCwsKwb98+scuySS16A57bqjaGplcfugPLH7yDoYmIiGyWzQWnjIwMTJ06Ff7+/pBIJNi9e3ebY1QqFUJCQuDs7IyoqChkZmaafP5Ro0bh4MGDqKqqwo0bN3Dw4EGUlpaa8RP0DI0teizZlIMv867AUSbBv36nwJMTBoldFhERUZfYXFedVqtFWFgYEhISMGvWrDavb926FUlJSUhJSUFUVBSSk5MRGxuLgoIC+Pj4AAAUCgV0Ol2b9+7fvx+hoaF49tlnce+998LT0xPjx4+HTNb+LNZNTU1oamoybtfU1JjhU9q2hmY9ntqYjUMXKiB3kOLD+eGYONxH7LKIiIi6zKbXqpNIJNi1axdmzJhh3BcVFQWlUok1a9YAAAwGA4KCgrB06VIsW7asw9d48sknMXPmTDz00EO3fP3111/HypUr2+zvqWvV1Ta24Im0bGT+VAlXJxnWxikRPbiv2GURERHdlqlr1dlcV93tNDc3IycnBzExMcZ9UqkUMTExOHbsmMnnuXbtGgCgoKAAmZmZiI2NbffY5cuXo7q62vhTXFzc+Q9g46rrW/D42kxk/lQJd7kDNj4RxdBERER2xea66m6noqICer0evr6+N+339fXFuXPnTD7P9OnTUV1dDTc3N6xbtw4ODu3/Z5LL5ZDL5VCpVFCpVNDr9Z2u35Zdr2vC/LWZyL9ag96ujtiYEMXZwImIyO7YVXAyl460Tv0sMTERiYmJxqa+nuR6XRMe+/gHFGhq4d1Ljk1PRmKEX8/rpiQiIvtnV1113t7ekMlk0Gg0N+3XaDTw8/Oz6LVVKhVCQ0OhVCoteh1rc0PbjHmftIYmH3c5ti4ez9BERER2y66Ck5OTE8LDw5Genm7cZzAYkJ6ejujoaIteOzExEfn5+cjKyrLodaxJ65imH3CurLWlactT4zG4Xy+xyyIiIrIYm+uqq6urQ2FhoXG7qKgIarUaXl5eCA4ORlJSEuLi4hAREYHIyEgkJydDq9UiPj5exKrtT3VDC+an/oAzV2rQ180JWxZFMTQREZHds7nglJ2djUmTJhm3k5KSAABxcXFIS0vDnDlzUF5ejhUrVqCsrAwKhQL79u1rM2Dc3HrS4PDaxhbEpWbiZEk1vNycsHnReAz1dRe7LCIiIouz6XmcrJGp80DYKm2TDnGpmci+dAO9XR2x+cnxCPW3v89JREQ9S4+cx4ksq7FFj0UbspF96QY8nB2w6YkohiYiIupRGJzMxN6fqtPpDXh2ywkc/fE63Jxk2PhEFEYF9KxpF4iIiNhVZ2b22FVnMAh4aUcePs8thZODFOvjIzkjOBER2RV21ZFZCIKAv3yVj89zSyGTSqB6bBxDExER9VgMTmZir111yd9eQNrRnwAAbz86BpNDLft0IhERkTVjV52Z2VNXXerhIvzlq3wAwMppIxH3mxBxCyIiIrIQdtVRl+w+UWoMTUmThzE0ERERgcGJbuHwhQq8vCMPABB/ZwiW3jtE5IqIiIisA4OTmdjLGKfTpdVYvDEbLXoBD4/pj9ceCoVEIhG7LCIiIqvAMU5mZstjnIor6zHrg6Mor21C9KC+SEtQQu4gE7ssIiIii+MYJ+qQSm0z4lIzUV7bhBF+7vhwQThDExER0S8wOBEamvV4Yn0WLlZoEdDbBWnxkfBwdhS7LCIiIqvD4NTD6Q0Clm7JxYnLVfB0ccT6BCX8PJ3FLouIiMgqMTj1cH/9Kh/fnr0GJwcp1sZFYIiPu9glERERWS0GJzOxxafqUg8XGWcF/+dvFYgI8RK3ICIiIivHp+rMzFaeqtt/pgyLN+VAEIBlU0bg6XsGi10SERGRaPhUHbXrVEk1nvtMDUEA5kYGYfHdg8QuiYiIyCYwOPUwJTfqkbA+Cw0tekwY6o2/TB/FCS6JiIhMxODUg9Q0tuCJtGzjXE3vzxsHRxn/FyAiIjIVvzV7CJ3egGc2n0CBphb93OVYu1AJd87VRERE1CEMTj3EG3vPIuN8OZwdpUiNUyKgt4vYJREREdkcBiczsebpCD794dJN0w6MDvQUtyAiIiIbxekIzMzapiM4WliBBamZ0BkEvHT/MDxz71CxSyIiIrI6nI6AUFShxZJPc6EzCJiu8EfipCFil0RERGTTGJzsVHV9C55Iy0J1QwsUQb3x90fGcNoBIiKiLmJwskM6vQGJm3NxsUILf09nfLQgHM6OMrHLIiIisnkMTnbojb1ncbiwAq5OMnwSp4SPu7PYJREREdkFBic7sy2r+P+eoJujQKi/+APUiYiI7AWDkx3JuVSJP+0+BQB4IWYYYkf6iVwRERGRfemxwWnmzJno06cPZs+e3ea1r776CsOHD8fQoUPxySefiFBdx12tbsDijblo0QuYMsoPS+/lE3RERETm1mOD03PPPYcNGza02a/T6ZCUlITvvvsOJ06cwFtvvYXr16+LUKHpGlv0eGpDDirqWtege/vRMEilfIKOiIjI3HpscJo4cSLc3d3b7M/MzMTIkSMREBCAXr16YcqUKdi/f78IFZpGEAQs23kSp0qr0cfVER8viICb3EHssoiIiOySVQanjIwMTJ06Ff7+/pBIJNi9e3ebY1QqFUJCQuDs7IyoqChkZmaa5dpXrlxBQECAcTsgIAClpaVmObclfJRxEbvVVyCTSvD+vHAEebmKXRIREZHdssqmCa1Wi7CwMCQkJGDWrFltXt+6dSuSkpKQkpKCqKgoJCcnIzY2FgUFBfDx8QEAKBQK6HS6Nu/dv38//P39zVZrU1MTmpqajNs1NTVmO/ev+f58Of627xwA4M9TQxE9uG+3XZuIiKgnssrgNGXKFEyZMqXd1999910sWrQI8fHxAICUlBTs3bsXqampWLZsGQBArVZ36tr+/v43tTCVlpYiMjKy3eNXr16NlStXdupaXXH5ej2e3XICggDMiQjC/PEDur0GIiKinsYqu+pup7m5GTk5OYiJiTHuk0qliImJwbFjx7p8/sjISJw+fRqlpaWoq6vD119/jdjY2HaPX758Oaqrq40/xcXFXa7h19Q36/DUxmzjcip/mTGSy6kQERF1A6tscbqdiooK6PV6+Pr63rTf19cX586dM/k8MTExyMvLg1arRWBgILZv347o6Gg4ODjgnXfewaRJk2AwGPCHP/wBffu23wUml8shl8uhUqmgUqmg1+s7/dlMIQgC/rjzFM6V1cK7lxwpj4dD7sDlVIiIiLqDzQUnc/n222/bfW3atGmYNm1ah86XmJiIxMRE1NTUwNPTs6vlteuTQ0X4Mu8KHKQSvD9vHPw8uZwKERFRd7G5rjpvb2/IZDJoNJqb9ms0Gvj5iTdTtkqlQmhoKJRKpcWucfhCBVZ/fRYA8NrDoYgc6GWxaxEREVFbNhecnJycEB4ejvT0dOM+g8GA9PR0REdHi1ZXYmIi8vPzkZWVZZHzF1fWY+mWXBgE4JFxgVgQzcHgRERE3c0qu+rq6upQWFho3C4qKoJarYaXlxeCg4ORlJSEuLg4REREIDIyEsnJydBqtcan7OxNs86AxRtzcKO+BaMDPLFq5igOBiciIhKBVQan7OxsTJo0ybidlJQEAIiLi0NaWhrmzJmD8vJyrFixAmVlZVAoFNi3b1+bAePdyZKDwx1lEswbH4x/p19AyvxwODtyMDgREZEYJIIgCGIXYU9+HhxeXV0NDw8Ps567oVkPFyeGJiIiInMz9fvb5sY49WQMTUREROJicDKT7niqjoiIiMTFrjozs2RXHREREVkGu+qIiIiIzIzByUzYVUdERGT/2FVnZuyqIyIisj3sqiMiIiIyMwYnIiIiIhMxOJkJxzgRERHZP45xMjOOcSIiIrI9HONEREREZGYMTkREREQmYnAiIiIiMhGDk5lwcDgREZH94+BwM6uurkbv3r1RXFzMweFEREQ2oqamBkFBQaiqqoKnp2e7xzl0Y009Qm1tLQAgKChI5EqIiIioo2pra28bnNjiZGYGgwFXrlyBu7s7JBKJcb9SqURWVlab403d/3MStpaWrPbq7u7zdeR9phz7a8fwPlrmfNZyH9t77Vb7rOle8j6a/po1/05ay33s6Hu7ei+t5T4KgoDa2lr4+/tDKm1/JBNbnMxMKpUiMDCwzX6ZTHbLm9nR/R4eHqL/cgPt19fd5+vI+0w59teO4X20zPms5T6299rtjreGe8n7aPpr1vw7aS33saPv7eq9tKb7eLuWpp9xcHg3SUxMNMt+a2Hu+jp7vo68z5Rjf+0Y3kfLnM9a7mN7r/E+mv99YtxHU68rFmu5jx19b1fvpa3dR3bV2QjOSG4feB/tB++lfeB9tA/deR/Z4mQj5HI5/vznP0Mul4tdCnUB76P94L20D7yP9qE77yNbnIiIiIhMxBYnIiIiIhMxOBERERGZiMGJiIiIyEQMTkREREQmYnAiIiIiMhGDkx0qLi7GxIkTERoaijFjxmD79u1il0SdNHPmTPTp0wezZ88WuxTqgK+++grDhw/H0KFD8cknn4hdDnUSf/9snyW+DzkdgR26evUqNBoNFAoFysrKEB4ejvPnz8PNzU3s0qiDDh48iNraWqxfvx47duwQuxwygU6nQ2hoKA4cOABPT0+Eh4fj6NGj6Nu3r9ilUQfx98/2WeL7kC1Odqh///5QKBQAAD8/P3h7e6OyslLcoqhTJk6cCHd3d7HLoA7IzMzEyJEjERAQgF69emHKlCnYv3+/2GVRJ/D3z/ZZ4vuQwUkEGRkZmDp1Kvz9/SGRSLB79+42x6hUKoSEhMDZ2RlRUVHIzMzs1LVycnKg1+sRFBTUxarpl7rzPlL36ep9vXLlCgICAozbAQEBKC0t7Y7S6X/w99M+mPM+muv7kMFJBFqtFmFhYVCpVLd8fevWrUhKSsKf//xn5ObmIiwsDLGxsbh27ZrxGIVCgVGjRrX5uXLlivGYyspKLFiwAB999JHFP1NP1F33kbqXOe4riY/30T6Y6z6a9ftQIFEBEHbt2nXTvsjISCExMdG4rdfrBX9/f2H16tUmn7exsVGYMGGCsGHDBnOVSrdhqfsoCIJw4MAB4ZFHHjFHmdRBnbmvR44cEWbMmGF8/bnnnhM+/fTTbqmXbq0rv5/8/bMenb2P5v4+ZIuTlWlubkZOTg5iYmKM+6RSKWJiYnDs2DGTziEIAhYuXIh7770X8+fPt1SpdBvmuI9kfUy5r5GRkTh9+jRKS0tRV1eHr7/+GrGxsWKVTLfA30/7YMp9tMT3IYOTlamoqIBer4evr+9N+319fVFWVmbSOY4cOYKtW7di9+7dUCgUUCgUOHXqlCXKpXaY4z4CQExMDB599FH85z//QWBgIP9RF5kp99XBwQHvvPMOJk2aBIVCgRdffJFP1FkZU38/+ftn3Uy5j5b4PnTo0rvJKt11110wGAxil0Fm8O2334pdAnXCtGnTMG3aNLHLoC7i75/ts8T3IVucrIy3tzdkMhk0Gs1N+zUaDfz8/ESqijqK99E+8b7aB95H+yDWfWRwsjJOTk4IDw9Henq6cZ/BYEB6ejqio6NFrIw6gvfRPvG+2gfeR/sg1n1kV50I6urqUFhYaNwuKiqCWq2Gl5cXgoODkZSUhLi4OERERCAyMhLJycnQarWIj48XsWr6Jd5H+8T7ah94H+2DVd5HszybRx1y4MABAUCbn7i4OOMx7733nhAcHCw4OTkJkZGRwvHjx8UrmG6J99E+8b7aB95H+2CN95Fr1RERERGZiGOciIiIiEzE4ERERERkIgYnIiIiIhMxOBERERGZiMGJiIiIyEQMTkREREQmYnAiIiIiMhGDExEREZGJGJyIiIiITMTgRERERGQiBiciIhHMnDkTffr0wezZs8UuhYg6gMGJiEgEzz33HDZs2CB2GUTUQQxORGS1Jk6ciOeff75Dx0skEkgkEqjVaovVZQ4TJ06Eu7v7LV976aWXMGPGjHbfu3DhQuPn3L17t2UKJKJbYnAiIruyaNEiXL16FaNGjTLuu+eeeyCRSPDmm2/edKwgCIiKioJEIsFf/vKX7i61XWq1GmPGjGn39X/961+4evVqN1ZERD9zELsAIiJzcnV1hZ+fn3FbEAScOHECAwYMwKlTp246dv369bhy5QoAYNy4cWarQaFQQKfTtdm/f/9++Pv7/+r78/LysGTJknZf9/T0hKenZ5dqJKLOYYsTEdmMvXv3wtPTE59++qnJ77lw4QJqa2sRFxd3U3Cqra3F8uXLsXDhQgBAeHg4AKCsrAwSiQT/+te/MHbsWDg7O2PkyJE4fPjwTee9fPky4uLi4OvrCxcXF4SFhRmPUavVOH36dJsfU0JTSUkJKioqAACTJ0+Gq6srhg8fjh9++MHkz0xElsPgREQ2YfPmzZg7dy4+/fRTzJs3z+T35eTkwNXVFXPnzkVBQQGam5sBAH/9618RERGBfv36wc/PD/379wcA49io1NRUJCcnQ61WIzg4GPPmzYPBYAAAXLp0CZGRkWhoaMAXX3yBkydP4plnnoGHh0eXP+fP11epVHjllVeQl5eH4OBgLFu2rMvnJqKuY1cdEVk9lUqFP/3pT/jyyy9xzz33dOi9ubm5GDNmDIYPHw5nZ2ecO3cOLi4u+OCDD5Cbm4tVq1bd1E2Xl5cHR0dH7NmzByEhIQCAN954AxERESgtLUVQUBCWLFmC8ePHY9u2bcb3DR06tEN1xcTEIC8vD1qtFoGBgdi+fTuio6OhVqvh5eWFbdu2wdvbGwAwbdo0fPjhhx06PxFZBoMTEVm1HTt24Nq1azhy5AiUSmWH35+bm4tx48ZBIpFgzJgxOHXqFLZs2YIlS5Zg6NChyMnJwcyZM43Hq9VqzJo1yxiaANzUknTp0iV8/fXXOHHiRJc+17fffnvL/Wq1GtOnTzeGJgAoKirCkCFDunQ9IjIPdtURkVUbO3Ys+vXrh9TUVAiC0OH3/xycgNZB28nJycjOzsZrr72GxsZGnDt37qYWJ7VaDYVCcdM5jh07Bm9vbwQEBECtVsPJyanNMeaiVqsxfvz4NvssdT0i6hgGJyKyaoMHD8aBAwewZ88eLF26tEPvvXjxIqqqqozBaOzYscjOzsbq1avh7u6OvLw86HQ648DwhoYGXLhwAXq93ngOg8GA5ORkxMXFQSqVwtHRETqdDvX19eb7kP9VW1uLixcvYuzYsTftZ3Aish4MTkRk9YYNG4YDBw5g586dHZoQMycnB05OTsY5neLi4lBeXm58ki43Nxf9+vVDUFAQAODUqVOQSCTYtGkTjh07hrNnz2LOnDmoqqrCq6++CgCIioqCp6cnlixZgrNnzyI/Px8pKSm4cOFClz9nXl4eZDIZRo8ebdx36dIl3Lhxg8GJyEowOBGRTRg+fDi+++47bNmyBS+++KJJ78nNzcWoUaPg6OgIAHB0dIS3tzckEonx9f9t3VGr1RgxYgReeeUVPPLII4iIiIBer8f333+P3r17AwD69u2LL7/8EhcuXIBSqcRdd92FL774Aj4+Pl3+jGq12jiI/WcnTpxA7969bxpzRUTikQidGTRARGSFJk6caBzH1BmJiYm4ceMGNm/ebN7CLEQikWDXrl23XZ6FiMyLLU5EZFfef/999OrVq80s4ab4taVOrMXTTz+NXr16iV0GUY/EFicishulpaVoaGgAAAQHB8PJycnk9wqCAE9PT3z22Wd48MEHLVWiWVy7dg01NTUAgP79+8PNzU3kioh6DgYnIiIiIhOxq46IiIjIRAxORERERCZicCIiIiIyEYMTERERkYkYnIiIiIhMxOBEREREZCIGJyIiIiITMTgRERERmYjBiYiIiMhEDE5EREREJmJwIiIiIjIRgxMRERGRif4/4sCziDWVA3AAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(hm.k_hm, hm.power_auto_tracer * hm.mean_tracer_den_unit**2 * hm.k_hm**3 / (2 * np.pi**2))\n", "plt.xscale(\"log\")\n", "plt.yscale(\"log\")\n", "\n", "plt.xlabel(\"k [$Mpc^{-1} h$]\")\n", "plt.ylabel(r\"$\\Delta^2(k) \\ \\ [{\\rm K}^2]$\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From v2.0.0, all of these HI component models, from [arxiv:2010.07985](https://arxiv.org/abs/2010.07985) are available in `halomod`. To use this model in full, simply call:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "hm = TracerHaloModel(\n", " hod_model=\"Spinelli19\",\n", " tracer_concentration_model=\"Maccio07\",\n", " tracer_profile_model=\"PowerLawWithExpCut\",\n", " z=1, # for default value of hod parameters,\n", " transfer_model=\"EH\",\n", ")" ] } ], "metadata": { "kernelspec": { "display_name": "cosmosis", "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.11.13" } }, "nbformat": 4, "nbformat_minor": 4 }