{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# kSZ Effect\n", "\n", "In this notebook we show how to use the `glasz` package to model the kSZ Effect profile as \n", "measured from [Schaan et al 2021](https://arxiv.org/abs/2009.05557). " ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "# preamble\n", "from __future__ import annotations\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pyccl as ccl\n", "from init_halo_model import ( # halo model\n", " a_arr,\n", " bM,\n", " cM_relation,\n", " cosmo,\n", " hmc,\n", " hmd,\n", " k_arr,\n", " r_arr,\n", ")\n", "\n", "import glasz\n", "\n", "# CMASS PARAMETERS\n", "z_lens = 0.55 # Mean z for CMASS\n", "a_sf = 1 / (1 + z_lens)\n", "\n", "# constituent fractions\n", "fb = cosmo[\"Omega_b\"] / cosmo[\"Omega_m\"] # Baryon fraction\n", "fc = cosmo[\"Omega_c\"] / cosmo[\"Omega_m\"] # CDM fraction\n", "\n", "# define 2-halo term\n", "xi_mm_2h = glasz.profiles.calc_xi_mm_2h(\n", " cosmo, hmd, cM_relation, hmc, k_arr, a_arr, r_arr, a_sf\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create Profiles\n", "\n", "We first create the gas (baryon) profile which we will then use to create the $T_{\\rm kSZ}$ signal. For more information on how to do this, visit the `profiles` notebook." ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "# Halo Mass\n", "M_halo = 3e13\n", "\n", "# 2-halo Amplitude\n", "A_2h = 1.0\n", "\n", "# GNFW Parameters\n", "alpha = 1.0\n", "beta = 3.0\n", "gamma = 0.2\n", "x_c = 0.5\n", "\n", "Rb = 10 * (hmd.get_radius(cosmo, M_halo, a_sf) / a_sf)\n", "\n", "# COMPUTE GNFW AMPLITUDE\n", "prof_nfw = ccl.halos.HaloProfileNFW(\n", " mass_def=hmd, concentration=cM_relation, truncated=False, fourier_analytic=True\n", ")\n", "\n", "prof_baryons = glasz.profiles.HaloProfileGNFW(\n", " hmd,\n", " rho0=1.0,\n", " alpha=alpha,\n", " beta=beta,\n", " gamma=gamma,\n", " x_c=x_c,\n", ")\n", "\n", "prof_baryons.normalize(cosmo, Rb, M_halo, a_sf, prof_nfw)\n", "\n", "\n", "# COMPUTE 3D DENSITY PROFILES\n", "def rho_2h(r):\n", " return (\n", " xi_mm_2h(r)\n", " * bM(cosmo, M_halo, a_sf)\n", " * ccl.rho_x(cosmo, a_sf, \"matter\", is_comoving=True)\n", " * A_2h\n", " )\n", "\n", "\n", "prof_baryons.rho_2h = rho_2h # add 2-halo term to baryon profile\n", "\n", "prof_matter = glasz.profiles.MatterProfile(\n", " mass_def=hmd, concentration=cM_relation, rho_2h=rho_2h\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From here we use the `kSZ` subpackage of `glasz` to compute the $T_{\\rm kSZ}$ signal." ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "# COMPUTE kSZ PROFILE\n", "theta = np.geomspace(0.5, 6.5, 20)\n", "\n", "\n", "def rho_gas_3D(r):\n", " \"\"\"\n", " We need to be careful about comoving units here. The kSZ code assumes\n", " that the density profile is in physical units. The array which it will\n", " be feeding into this function is in physical units r = aχ. The CCL\n", " profile assumes comoving units so we need to convert the input to\n", " comoving units before passing it to the profile. Then we need to convert\n", " the density profile back into physical units before returning it by\n", " dividing by the scale factor a^3.\n", " \"\"\"\n", " return (fb * prof_baryons.real(cosmo, r / a_sf, M_halo, a_sf)) / a_sf**3\n", "\n", "\n", "T_kSZ_150 = glasz.kSZ.create_T_kSZ_profile(theta, z_lens, rho_gas_3D, \"f150\", cosmo)\n", "T_kSZ_090 = glasz.kSZ.create_T_kSZ_profile(theta, z_lens, rho_gas_3D, \"f090\", cosmo)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can plot the difference between the 2 frequency bands for a given profile below." ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAisAAAF1CAYAAAApwqoyAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVepJREFUeJzt3Xl4lNX9/vH3ZN8D2ROSsMu+gwi4sCMqyiJirRT3WqlKsZW60K9ahKKtYiuiaFXUqvyqgloXRBRQXNj3HYIBkhBIyL7OzPP748kKAZKZCZkk9+u6cmXmmZkzJ0GY23M+5xyLYRgGIiIiIm7Ko6E7ICIiInI+CisiIiLi1hRWRERExK0prIiIiIhbU1gRERERt6awIiIiIm5NYUVERETcmldDd6Axs9vtpKSkEBwcjMViaejuiIiINBqGYZCbm0tcXBweHucfO1FYcUJKSgoJCQkN3Q0REZFG6+jRo8THx5/3OQorTggODgbMX3RISEgD90ZERKTxyMnJISEhoeKz9HwUVpxQPvUTEhKisCIiIuKA2pRRqMBWRERE3JrCioiIiLg1hRUHLFy4kK5duzJgwICG7oqIiEiTZzEMw2joTjRWOTk5hIaGkp2dfd6aFZvNRmlp6UXsmTjC29sbT0/Phu6GiEizUNvPUFCBbb0yDIO0tDSysrIauitSSy1atCAmJkb75oiIuBGFlXpUHlSioqIICAjQB6AbMwyDgoIC0tPTAYiNjW3gHomISDmFlXpis9kqgkp4eHhDd0dqwd/fH4D09HSioqI0JSQi4iZUYFtPymtUAgICGrgnUhflf16qMRIRcR8KK/VMUz+Ni/68RETcj8KKiIiIuDXVrIiIiEi9MgzILISUPEjJNb8Op9X+9QorIiIi4pSC0soQUu2rSjgptlV/jb2o9u1rGkiqWbt2LePGjSMuLg6LxcLy5cvPes5tt92GxWKp9nXZZZdVe05xcTH3338/ERERBAYGcv3113Ps2LELvn9aWhoPPvggHTp0wM/Pj+joaC6//HJefvllCgoKKp7Xpk0bFixYcNbrn3jiCXr37l3XH1tERM7BaofjObDhOHy8DxZthMe/hTs/gbH/gV6vQJeXYMTbMHU5zFoFL6yH/+6BdUchKasyqEQFQu9ouKYDTO1Z+z5oZEWqyc/Pp1evXtx+++1MmjTpnM+7+uqreeONNyru+/j4VHt8xowZfPrpp7z//vuEh4fz0EMPcd1117Fp06ZzLgk+fPgwQ4YMoUWLFsydO5cePXpgtVrZv38/r7/+OnFxcVx//fWu+UFFRATDgNNFcDwXUnPLvpeNhpRfO5EP9lrsdR/iA7HBEBcMcUFl36t8RQcaFBfbOZFpIy3DStLRnFr3U2HlIjIMKLRe/Pf194LaLnIZO3YsY8eOveDzfH19iYmJqfGx7Oxs/v3vf/P2228zcuRIAN555x0SEhL4+uuvGTNmTI2vu++++/Dy8mLjxo0EBgZWXO/RoweTJk3CkZMhalrd07p1a44cOVLntkREGpvy6ZmqQeR4buW1lDwoqsXnkrdHWRCpIYTEBZmPBfsYnM61cyLTWhFIjiZZ2Zhh3j+RaaWgqPLfcWtJbq1/DoWVi6jQag6VXWx77oMAb9e2uXr1aqKiomjRogVXXXUVTz/9NFFRUQBs2rSJ0tJSRo8eXfH8uLg4unfvzg8//FBjWMnIyOCrr75i7ty51YJKVY4sK05NTa24nZ+fz9VXX82gQYPq3I6IiLux2uFEXmVdSHkgqVorklXLupDIgOrh48xAEhEAGAYZOTZOZJjB48RxK1sybXyZaSUtw0p6po3i0gv/T2XLEA+iw7wI9Q/gu1r+rAorUmdjx45l8uTJtG7dmqSkJGbPns3w4cPZtGkTvr6+pKWl4ePjQ8uWLau9Ljo6mrS0msu/Dx48iGEYdOrUqdr1iIgIiorMv23Tp09n/vz5FY/NmjWLxx9/vNrzS0pK6Nq1a8X98tEfwzCYNGkSoaGhvPLKK47/8CIiF0H59ExNharlX7WdngnyOXcIaRUM0YHg6wU2m8HJLBsnMqykZVo5ccDGrvLbmTbSM61Ybed/L4sFwkM9iQn3IjrMk5gwL6LLbkeHexHd0hNfH7NcNicnh789ULvfh8LKReTvZY5yNMT7utKUKVMqbnfv3p3+/fvTunVrPvvsMyZOnHjO1xmGccHRkTMfX79+PXa7nV//+tcUFxdXe+xPf/oTt912W7Vr//znP1m7du1Z7T766KP8+OOPbNiwoWJbfRGRhlJiMwPHsZzK0ZDyGpGUOkzPeHlAbND5R0VCfMves9Tg5GkraZlmINlz0MrqsimatEwrp7Js2O3nfz8PD4hq6Ul0mFdFIImuEkiiWnrh7eX6zTUVVi4ii8X10zHuIDY2ltatW3PgwAHAHM0oKSnh9OnT1UZX0tPTGTx4cI1tdOjQAYvFwt69e6tdb9euHUCNASMiIoIOHTpUuxYWFnbW89555x2ef/55Vq9eTXx8fN1+OBERB5TvK5KcA0ezITnbvJ2cbd5PyavdqEhkwBlh5IyvyADwKMsGRSX2iimatKNWvisLISfKAklmjo0Llf55e0FUyypBJNzLDCNh5mhJRKgnnp4Xf6dvhRVxWkZGBkePHq04qbhfv354e3uzcuVKbrrpJsCsHdm5cyfPPPNMjW2Eh4czatQoXnzxRe6///5z1q3U1Y8//shdd93FK6+8ctbyahERZxRZzZGRo2UhJDm7+u38Cxwx5usJCaHmVExNIyMxQeBX5VM6v9BOWoaVE5lWjhy38VPZ7fKAkpV3gWERwNfbQkx49dGQmCqBJCzEEw8P9zt2RGFFqsnLy+PgwYMV95OSkti6dSthYWEkJiaSl5fHE088waRJk4iNjeXIkSM8+uijREREMGHCBABCQ0O58847eeihhwgPDycsLIw//vGP9OjRo2J1UE1eeuklhgwZQv/+/XniiSfo2bMnHh4ebNiwgb1799KvX786/SxpaWlMmDCBm2++mTFjxlTUy3h6ehIZGenAb0dEmhPDgPSCyhBSMUJSNkqSlnfhNmKCIDHEDCWJIZAYWnY7tPqoiGEY5OSXLes9ZWXD/rIVNGWBJC3DSl7hhYdiAv0sNQaR8oASGuTRKM9AU1iRajZu3MiwYcMq7s+cOROAadOm8eabb+Lp6cmOHTt46623yMrKIjY2lmHDhrF06VKCg4MrXvf888/j5eXFTTfdRGFhISNGjKh4/bm0b9+eLVu2MHfuXB555BGOHTuGr68vXbt25Y9//CP33Ve3gp+9e/dy4sQJlixZwpIlSyqua+myiJQrKDVDSMWIyBmjJBeqGwn0NoNHYijEl4WR8lASH1J9ZMRmM0jLsHL8pJUf9ls5nl5KyikraWUjI4XFFw4jIYEeFXUiMWcFEi+CAprmXq8Ww5HNKwQwK5lDQ0PJzs4mJCSk2mNFRUUkJSXRtm1b/Pz8GqiHUlf6cxNpWuyGOQJy5hRNeS3JyYLzv97DYk7NlI+GJIRUhpPEEAjzr76PlbVKIDmebuX4ydKK22kZVmwXmKkpX9YbE3bGFE24ubLG36/phJHzfYaeSSMrIiLSqOUWVylerToykg3Hcs2VN+cT4gOJLapM01QJJHHB4HPGgHCp1QwkBw5VCSNlX2kZ1vOuqPHxttAq0ou4SC/iy77HhHudtaxXqlNYERERt2a1m8t5q9WOVAklpy+w8ZmXh1nEmnhm3UjZ7dAaBlFLrQapGVY2HS0LJOmVgeRE5vkDiW/VQBLlRVykN/FRXrSK9CI81D0LWN2dwoqIiDQ4wzCnZA6dhsOnK78fPm2uuLFdoGAhzL/mIJIQai799aphwKKk1CD1lJVdVaZqzEBSSnqm7bxLi/18LNXCSKvykZIoc3lvYyxidWcKKyIictEUWSHpNBzOqh5IDp+G3JJzv87X0yxYrVYzElq50ibIp+bXlZQaHDtRfXQk5aSVY+mlpJ8+/74jfr7mCIn55V0WTCpHSBRILh6FFRERcSmjrKi1IoxkVY6WHM+Bc+UDD4sZSNq1gPYtoV1LaB8GbUIhOqhyme+ZikvsJKWUh5Cy72WjJScvEEj8fS20iioLI2UjI+YoiTctQxrnMt+mqNmHlQkTJrB69WpGjBjBBx980NDdERFpNApKq4+MHCr7SsoyHzuXEF8zjJQHknZlt1uHmmfU1KSoxE7qKTOMVB0dSTlp5WTW+QNJgJ+lYnSkVXkYifImLtKLlsEKJI1Bsw8rDzzwAHfccUe1fThERMRkN8zi1qphpDycpJ5nUzRPixk+2p0RSNq1hPAzlvtWlZ1n40hqKb+klprf00o5dsIMJOcT6GehVVRlGKkaTlo00o3QpFKzDyvDhg1j9erVDd0NEZEGlVdScyBJyjr/xmhh/ua0TbszRkoSQ89e8ltVdp6NX9JKOZJSGUqOpJZyOufcy2wC/S3EVxkdaRXlXRFMGuvOrFI7jTqsrF27lmeffZZNmzaRmprKsmXLGD9+fLXnvPTSSzz77LOkpqbSrVs3FixYwBVXXNEwHRYRaUA2u3my78HMKlM3Web39Pxzv87bA1q3KAsjLaqPlLS8wCHmjoSS6DBP2sR60ybWm9ax3iTGmKEkJFCBpLlq1GElPz+fXr16cfvttzNp0qSzHl+6dCkzZsyoOHPmlVdeYezYsezevZvExMQ6v19xcTHFxcUV93NycpzqvzvKzc1l9uzZLFu2jPT0dPr06cMLL7zAgAEDKp5jGAZPPvkkixcv5vTp0wwcOJCFCxfSrVu387adk5PDs88+y0cffcThw4cJCAigXbt2TJ48mbvvvrvihOahQ4fSu3dvFixYUO31b775JjNmzCArK8vVP7ZIk1JiMwPI3lOwP6NspCQLfsmC4vPMpkQGnD1t076lWfRa09LfqnLyzembIymVgaSuoaRNrDetY7yb1C6t4hqNOqyMHTuWsWPHnvPx5557jjvvvJO77roLgAULFrBixQoWLVrEvHnz6vx+8+bN48knn3S4v43BXXfdxc6dO3n77beJi4vjnXfeYeTIkezevZtWrVoB8Mwzz/Dcc8/x5ptvcskllzBnzhxGjRrFvn37qp0PVFVmZiaXX345OTk5/PWvf6Vfv374+Phw8OBB3n33Xd59912mT59+MX9UkUbPMOBEvhlK9p6CPadgX4Y5clJ6jozg6wltWlQPI+1aQtuWEOp74fcsDyXlNSXltzNrEUrKA4lCidRVow4r51NSUsKmTZv485//XO366NGj+eGHHxxq85FHHqk42A/MkYKEhASn+ulOCgsL+fDDD/n444+58sorAXjiiSdYvnw5ixYtYs6cORiGwYIFC3jssceYOHEiAEuWLCE6Opp3332X3/72tzW2/eijj5KcnMy+ffsqQg9A586due6663DkiKo2bdrwyy+/nHVdx11JU1RQao6SlAeSPWUBJescu7cG+0CncOgcUb2WpFUweNYiI7gqlCTGeBOgUCJOarJh5dSpU9hsNqKjo6tdj46OJi0treL+mDFj2Lx5M/n5+cTHx7Ns2bJqUx5V+fr64utbi//1OAfDMCgqufgfpH4+llrN81qtVmw221kH+Pn7+/P9998DkJSURFpaGqNHj6543NfXl6uuuooffvihxrBit9tZunQpt956a7WgUpUj89AbNmzAZjPHtG02GzfeeCPe3t51bkfEndgNcwv5Padg36nKUPJLds37k3hazFGRLhHQuSycdI4wQ0lt/lrl5NvOCiRHahFKqo2SlI2UKJRIfWmyYaXcmR+ChmFUu7ZixYqL1peiEoNr/3Dsor1fuc+ej8ff98L/agUHBzNo0CD++te/0qVLF6Kjo3nvvff4+eef6dixI0BF0KspBNY0ygFw8uRJsrKy6NSpU7Xr/fr1Y9++fQCMGzeO9957r+Kxl156iddee63a861Wa7UgFRkZWXH7wQcfJDU1lQ0bNlzw5xRxF1lFldM3e8tGTPZlnHuPksgAM4h0Ci8LJxHQIQz8avEveW6BnSMpJWeNliiUSGPQZMNKREQEnp6e1UZRANLT08/6oJVKb7/9NnfccQetWrXC09OTvn37csstt7B58+Zqz7tQCKzJmY8vW7aMkpISZs2aRWFhYbXHfv3rX/PYY49Vu/bRRx8xd+7cs9pdvHgx//73v1m3bl21ACPiLsoLXstDyd5TsDfD3OW1Jr6e0DEculQZKekcAREBF34vu93g2EkrB5JL2J9cwqFjJQol0ug12bDi4+NDv379WLlyJRMmTKi4vnLlSm644Qan2l64cCELFy6smIKoLT8fC589H+/UezvCz6f2Uyzt27dnzZo15Ofnk5OTQ2xsLFOmTKFt27YAxMTEAOYIS2xsbMXrzhcCIyMjadGiBXv37q12vXxFVnBw8FkrfEJDQ+nQoUO1a1FRUWe1vXr1au6//37ee+89evXqVeufU6Q+lBe8Vgslp8zVOOcqeI0PqRwlKZ/GadPiwqtvAGx289yb/ckl7D9awoHkEg4eK6GgqObp5qhzrL5RKBF316jDSl5eHgcPHqy4n5SUxNatWwkLCyMxMZGZM2cydepU+vfvz6BBg1i8eDHJycnce++9Tr3v9OnTmT59Ojk5OYSGhtb6dRaLpVbTMe4gMDCQwMBATp8+zYoVK3jmmWcAaNu2LTExMaxcuZI+ffoAZjHzmjVrmD9/fo1teXh4cNNNN/HOO+8we/bsc9at1NXBgweZNGkSjz76aEWxr8jFUlBqTtlUDSV7TkF2cc3PD/Y5ewqnUzgE17IMzmYzSD5RWjFisv9oKQePlVBUfHYw8fW20D7em46JPlyS4EObOIUSadwadVjZuHEjw4YNq7hfvlJn2rRpvPnmm0yZMoWMjAyeeuopUlNT6d69O59//jmtW7duqC67vRUrVmAYBp06deLgwYP86U9/olOnTtx+++2AGbhmzJjB3Llz6dixIx07dmTu3LkEBARwyy23nLPduXPnsnr1agYOHMhTTz1F//79CQwMZPv27fz444907969Tv0sLCxk3Lhx9O7dm3vuuafadF/56I+IqxRZYfdJ2HYCtqaZ349k1a7gtUsEdKpDwSuYweSXtFIzlFRM55RSXHr2O/r5WOiQ4MMlCWXhJNGHxGhvPD0bx/8YidRGow4rQ4cOveAy1fvuu4/77rvvIvWo8cvOzuaRRx7h2LFjhIWFMWnSJJ5++ulqq2wefvhhCgsLue+++yo2hfvqq6/OuccKQHh4OOvXr2f+/Pk8++yzJCUl4eHhQceOHZkyZQozZsyoUz9PnDjB3r172bt3L3FxcdUe09JlcYbNbk7bVA0me06BtYZpnPKC16pTOLUteC1ntRkcSTGDyYGjZcHkeCklNQQTf9+yYFI2YtIx0YeEaC88z3UcsUgTYTH0L7vDyqeBsrOzCQkJqfZYUVERSUlJtG3b9qylwOK+9OfWvBiGWeS69QRsSzO/70g3z8k5U7g/9IqG3jHm9+5RtSt4rarUapBUHkzK6kwOHy+htIazdwL8LHQsDyaJZjCJj/TCQ8FEmojzfYaeqVGPrDQURwtsRaRhZRfD9rJQsv2E+b2mM3H8vaBndGU46RkN8XWYxgEoKTVISimpNpWTlFKKtYZ/NgL9LWYgqRJO4iIUTETKKaw4wNECWxG5eIqs5vRN+VTOtjTzfJwzeVrMmpLeVcJJh7DarcYpV1xi5/Dx0ooVOeXBxFbD1FFwgEdZMPGuGDGJi/DSAX0i56GwIiKNnt0oqzNJq5zS2XOq5uXCiaGVwaRXDHSPBP86bHxstRkcOFrC3iOVNSZHUkux1/BeIYEeldM4ZaMmMeGeCiYidaSwIiKNTlpe5YjJ1jSzziS3hjqTsPI6k7Jg0ivavFYXhcV29iSVsP1gETsOFbMnqaTGYzNaBHlUjJSUh5PoMAUTEVdQWBERt5ZTbNaXbDtRGU5OnKPOpEdUZSjpFQ0JIXWrMwHIzrOx81Ax2w8Ws+NQMQeSS86azgkJ9KBLm8r6kksSfYhooWAiUl8UVuqZvaaxYXFb+vNqeDnFsOE4/HgcfjoGu06a0zxVeVjMpcLlwaR3tLk9fV3qTMqdyLSy42AxOw4Ws/1QMb+knn0wT1RLT3p08KVHe196dPCldYy3il9FLiKFFQfUZjWQj48PHh4epKSkEBkZiY+Pj/6vy40ZhkFJSQknT57Ew8MDHx+fhu5Ss5FbDBtS4MdjZjjZWUM4iQ+BPlWCSbcoCHDggG3DMPglzcqOg0UV4SQ98+y/x61jvOjRwa8inMSE659KkYakfVaccKE14iUlJaSmplJQUNAAvRNHBAQEEBsbq7BSj/JKzHDy0zEzoOxIPzuctAmFQQlwWSu4LB5ighx7r/Ji2PKRkx2HisnJrz565uEBHRN86Fll5CQ0yNPBn05Eakv7rLgJHx8fEhMTsVqt2pOlEfD09MTLS0tIXS2/Sjj56bhZf2I7I5y0DjVDyaB4M6DEnnsz5PMqKrGzO6k8nBSx+8jZZ+f4elvo0rYsnHTwo2sbH/x1Zo6IW1NYqWcWiwVvb+9q29WLNGX5JbAp1Rw1+fFYzeEkMdQMJYPizZAS52A4ycm3seNQccXIyf4aimGDAzzo1s6Hnh386NnRl44JPnh7KZCKNCYKKyLilIJS2JRiFsSWh5Mzz9GJDzGDyaB4GNjKvO+Ik6et5iqdsimdpJSzi2EjWpjFsD3LpnTaxKoYVqSxU1gRkTopLK0cOfnpmLmc+MzN11oFVwkn8eYSYkeczLKyYVcR2w6Ye5ykZZw9nZoQ7VUlnPhp0zWRJkhhRUTOq8hqhpPygtitaWeHk7igyoLYQQmOhxObzWBXUjE/7yzi512FHD5efeTEwwIdEnyqLSNuGaxiWJGmTmHFATrIUJq64znwdRJ8fdgsii054z/12KDKepNB8Y5tvlYuM8fGht2F/LyziI17CskrrCxwsVigSxsf+nY2lxF3a+dLgIphRZodLV12Ql2WXYm4M7sBO9Nh5WEzpOw+Wf3x6MDKaZ1B8WaBrKPhxGY32PdLCT/vLOTnXUXsT66+T35IoAcDuvoxsJs/A7r6aRmxSBOlpcsickFFVlh3FFaVBZSqW9h7WKBfLIxsByPbQvuWjocTMLew37iniJ93FrJ+d9FZe51ckujDwG5+DOzuT6fWPniqIFZEqlBYEWlGTubDN0dgVRKs/QUKrZWPBXrDla1hVDsY1qbuB/5VZbcbHDxWWjZ6UsieIyVUHcMN9LfQv4s/A7v5cWlXf8JCNXoiIuemsCLShBkGHMisnN7ZkgpV531jg8zRk1HtzOJYXyf+RcgrsLNxb/noSSGnc6qPnrRr5c3Abv4M7O5H17a+eHlq9EREakdhRaSJKbXB+hRzemdlEiRnV3+8R5QZTka0hW6Rjk/vGIbB4eOl/LzLrD3ZdbiYqudA+vta6NfZj0u7mSMokS31z42IOEb/eog0AdnFsOaIOYKy+gjkVKlZ9fWEwQmVAcXRc3YACorsbN5bVBFQTmVVXybUOsbLDCfd/enR3lc7xYqISyisiDRSydnm0uKVh82RlKq7xob7w/C2ZnHsFYkQ6OC5jIZhkJxm5eddhazfVcj2g8VYq+QTX28LfTr5MrC7PwO7+et0YhGpF/qXxQHaZ0UaglG2vPiLg+b0zv6M6o93DDPDych20CcGPJ3YjiT1lJVVG/JZtbGAX1Krb8zWKtKrYuVOr45++Hhr9ERE6pf2WXGC9lmRi+FEHizbCx/sMYtly3laYEArGFUWUNq0cO59TufaWLO5gK/X57M7qXIeydsLel9i7ntyaTc/4qN0KKeIOE/7rIg0coWlsOIQfLgHvj9qbtoGZv3JyHYwumx5caifk+9TZOf7bYWs2pDPxr1FFQWyHhbo08mP4QMCuKJ3AEH+2jVWRBqOwoqImzAMs/bkwz3w2QHIq1Ik2z8WbuwK13SEUF/n3qfUarBhdyGrNhbww7ZCiksrB1c7Jfow4tIAhvULJFx7n4iIm1BYEWlgydlmQPlwDxzNqbweHwKTOsOkLtC6hXPvYbcb7DxczKr1BazZUlBtB9lWkV6MGBDAiAGBJERrikdE3I/CikgDyCk2R08+3AMbUiqvB/nANR3MgHJpK3M6xhmHjpWwamMB32zMJz2zsiA8LMSDYf0DGdE/gE6tfbA4s5e+iEg9U1gRuUhsdvgu2QwoKw5BcVl2sACXJ5oB5er24O/k4EZahpVvylbyJKVUruQJ8LNwRe8ARl4aSO9LfHX+jog0GgorIvVsf4a5kmfZXkivclhghzC4sQtM6OzcRm1gHhS4elMBqzYWsPNQccV1by8Y2M2fEQMCuay7H74+KpQVkcZHYUWkHmQUwCf7zVGUHemV11v4wQ2dzFqUntHOnWRcWGznh+2FfL0+n417irCVlaFYLNCroy8jBwRyRZ8AggMUUESkcVNYEXGRUpt5mvGHe8yTjct3lPXygOFtzGme4W3Bx4lFNlabwcbdRazamM+6bYUUlVSu5OmY4M2IAYEM6x9AZAv91RaRpkP/ojlAO9hKVbnF8N4ueH0LpOZVXu8RZQaU6y+B8ADn3iM7z8b/vs9j+Zo8MrIr/7uLjfBi5IAAhvcPpHWsVvKISNOkHWydoB1sm7eUXHhjK7y3E3LL9kSJ8IdJXc1alEvCnX+P5LRSPvwml69+zq/YD6VFkAfD+gcwckAgndtoJY+INE7awVakHu0+CYs3w6f7K6d6OoTB3X1gfGfwc/JvlWEYbNpbxAff5LJ+V1HF9Q7x3tw4PJih/QJ1Ho+INCsKKyK1YBjmsuPFm83v5S6Lh3v6mlvfO7sSuKTU4Ov1+Xz4bW7FkmOLBQb38OfG4cH07OirURQRaZYUVkTOo8QG/9tvhpQ9p8xrHha4tqMZUnpGO/8emTk2Plmbyydr88jKM4dq/HwtjB0UyMShwbTSwYEi0swprIjUIKcY3t1p1qSklRXNBnjDzd3gjj6Q4IISpUPHSvjgm1y+2ZhPqdW8FhXmyYSrgrl2SBBBWnIsIgIorIhUk5ILr281i2bLDxKMDIA7esOvezh/yrHdbvDzriI++CaHLfsqN2/r2taHG4cHc0XvADw9NdUjIlKVwooIsOskLN4E/ztQWTTbMcyc6rmhE/g6+TelsNjOVz+Z9SjH0s1hFA8PuLJ3ADeOCKZrWyePUhYRacIUVqTZMgxYm2yGlO+PVl4fHA/39IOhrZ3bYRbg5Gkry9bk8dn3eeQWmCko0N/CtUOCmDA0mOgw/RUUEbkQ/UspzY7NDh/vg1c2wd4M85pnlaLZHi4omt1zpJgPvsllzeYC7GUjNa0ivZg4LJirLwvE30/1KCIitaWwIs3K2l/g6e8qQ0qgN9zc3axJiXeyaNZmN1i3rZAPvsmtdphg746+TBoRzGXd/XXSsYiIAxRWpFnYdwqe/h7W/GLeD/WFe/vBr3uat521cU8hiz7MqtgfxcsThvcPZNLwYDom+Dj/BiIizZjCijRp6fnw3E+wdBfYDfD2gN/0ggcuNU9AdtYvqaW8/NFpfi7baTY4wIMbrgrihiuDCQ914sRCERGpoLDiAB1k6P4KS+HVLfDyRsg3BzsY2wH+PATatHC+/ew8G29+ls2n3+Vht4OnB4wfGszUsSGEBCqkiIi4kg4ydIIOMnQ/dgM+2gPP/li5mVufGHjsChgQ53z7JaUGy9fk8vYX2eQXmn91hvT0554JLUiI1k6zIiK1Ve8HGX7yySd1fs2oUaPw9/d35O1EauWHozDnO3PPFID4YJg1BMZd4vwSZMMw+G5rIYuXZ5Fy0twnpX28N7+b1JK+nVwwnyQiIufkUFgZP358nZ5vsVg4cOAA7dq1c+TtRM7rYCbM+x6+TjLvB/vA7wfAbb2dPwEZYH9yCS99cJrtB80VPmEhHtx5fQtGXxao1T0iIheBw/+Up6WlERUVVavnBgcHO/o2IueUUQDP/wzv7gCbYe6VcmtPmDEQwlwwiHcyy8rrn2Tz1c/5GAb4eFuYMjKYm0eFaJ8UEZGLyKGwMm3atDpN6dx6662q6RCXKbKaBwwu3AC5Zef3jGoHj1wO7Vs6335hsZ2lK3P4f1/nUlRi1qWMHBDAXTe0IEo7zoqIXHQqsHWCCmwvLrsBn+yDZ36A47nmte5R8PgVMCjeBe3bDVauz+e1j7PJyDZXenVv78vvJrWgSxud3SMi4kr1XmArcrEdzIQ/roQtaeb92CB4eDCM7wyuKBvZdqCIlz44zYGj5jrnmHBP7pnQkqv6+GNxtjpXREScUuewUlhYSGZmJq1atap2fdeuXXTr1s1lHRMBczTlre0w9zsotpnb4/+uP9zVB/xdsFL4+MlSFi/L4ruthQAE+Fm49epQJg4LxsdbIUVExB3UKax88MEH/OEPfyAsLAzDMHj11VcZOHAgAFOnTmXz5s310klpntLy4E8rzZORAa5MhGdHQUyQ823nFdh5+4tslq3OxWozR2euvTyI264LpWWwNnUTEXEndQorc+bMYfPmzURGRrJx40amTZvGY489xi233IJKX8SV/rcfHv0GsovB19Pc1O03PZ3fLwVg3fYC/v5OJtl55nHIA7r6ce/EFrSN0xk+IiLuqE5hpbS0lMjISAD69+/P2rVrmThxIgcPHtS8vrhEdjH85VtYvs+83zMKnh8DHcKcb7uk1ODlj06zfI25tW3rGC9+N6kll3bTZoUiIu6sTmElKiqK7du307NnTwDCw8NZuXIl06ZNY/v27fXSQWk+1h2Fh76C1DxzWub3A8wDB71dMCtzJLWUOf8+xeGyU5EnjwjmzutbqC5FRKQRqNPS5WPHjuHl5UVMTMxZj61bt44hQ4a4tHPuTkuXXaPIai5H/vcW836bUHM0pW+s820bhsFn6/JZ+N/TFJcatAz2YNZvwjWaIiLSwOpt6XJ8/Lk3s2huQUVcY2c6zFgBBzLN+7f2MOtTAlyw0ie3wM4//pPB2i3mSp9+nf14ZFo4YaEqoBURaUxcss/KqlWrWLVqFenp6djt9mqPvf766654C2libHZ4eRM8/xOU2iEyAOaPhBFtXdP+zkPFzHnjFOmZNjw94M4bWnDTiGA8dJaPiEij43RYefLJJ3nqqafo378/sbGxKrSVC0rOhj+sgI2p5v0x7WHecAgPcL5tm93g3S9zWPJZNnYD4iK9ePz2cDprB1oRkUbL6bDy8ssv8+abbzJ16lRX9EeaMMOApbvgqbWQXwpBPvDkVTCpi2uWJJ88bWXumxlsO2Cejjzy0gAenBJGoL8OHRQRacycDislJSUMHjzYFX1pNBYuXMjChQux2WwN3ZVGI7MQZn0NXx02718aB8+NgQQX1SV/v7WAv/8nk5x8O/6+Fh68OYzRAwNd07iIiDQopw8ynDVrFkFBQcyePdtVfWo0tBqodpJOw28+Nqd/vD3gj4Ph7j7g6YIBj+ISOy9/lMXHa829Uy5J9OHxO8KJj3JBha6IiNSbi3qQYVFREYsXL+brr7+mZ8+eeHtX/5B47rnnnH0LacQ2pcKdn8DpInMUZfF10DXSNW0npZQw5/UMksr2TrlppLl3ireX6qZERJoSp8PK9u3b6d27NwA7d+6s9piKbZu3Lw/CA1+aBxD2jILXr4dIF8zM1LR3yp+nhTOgq/ZOERFpipwOK99++60r+iFNzBtb4ck1YGAuR35xbP3snTKgqx+zfhNOWIj2ThERaapcss+KSDm7AXO+q9yN9tYe8ORQ8HJBfcqOg0U8/UYG6adteHnCXTe04Mbh2jtFRKSpcyiszJw5k7/+9a8EBgYyc+bM8z5XNSvNR5HV3D/l84Pm/T8PgXv7uWZZ8qoN+fxtSQY2O7SK9OLxO8Lp1Fp7p4iINAcOhZUtW7ZQWlpacftcVLPSfJwuhLs+NTd68/GEv4+CGzq5pu1P1ubywtLTGAZc1TeAP90aRoCf9k4REWkunF663Jxp6bIpORumLYfDWRDiC69eB5ed+xipWjMMg/dW5PDaJ9kA3HBlEPff1FLTPiIiTcBFXboM5vLl7du3n3U2kMViYdy4ca54C3FT29Lgjk/gVCG0CoY3b4BLwp1v1zAMFi/LYunXuQD8+uoQ7hgXqtE6EZFmyOmw8uWXXzJ16lQyMjLOesxisWiX1yZs1WGY/gUUWs29U968HqKDnG/XZjdY8F4mn63LB+DeiS24aWTzHbkSEWnunJ74//3vf89NN91Eamoqdru92peCStP1zna4639mULmqNfz3RtcElVKrwZzXM/hsXT4eFvjTrWEKKiIizZzTIyvp6enMnDmT6OhoV/RH3JzdgGfWwaJN5v0p3eDpYeDtgm1OCovtPPHqKTbsLsLLEx67PYKr+rrgKGYREWnUnB5ZufHGG1m9erULuiLuzmaHmV9VBpWZl8H8Ea4JKnkFdma9eJINu4vw87Hw9O8iFVRERARwwWqggoICJk+eTGRkJD169DjrbKAHHnjAqQ66s+a0Gsgw4P/WwJJt5gZv80fAjV1d03Zmjo1ZL6Zz6FgpQf4W5t4XRff22kNFRKQpu6irgd59911WrFiBv78/q1evrrZaw2KxNOmw0py8sskMKgAvjIHrLnFNu2kZVh7+VzrH0q20DPHgmd9H0T7exzWNi4hIk+B0WHn88cd56qmn+POf/4yHhzbqaoo+2gvz1pm3/3Kl64JKclopf/pnOiezbESHefL3B6JoFeWCA4RERKRJcTqslJSUMGXKFAWVJuq7X+BPK83bd/eFO/u4pt39ySXMejGd7Dw7rWO8eOaBKCJb6KgqERE5m9MJY9q0aSxdutQVfRE3szMdfvsZWO1w/SXw6OWuaXfbgSJmLjhBdp6dTok+LJgZraAiIiLn5PQnhM1m45lnnmHFihX07NnzrAJbHWTYOB3Ngds+hvxSGBRvnvXjil3uf9pRyBOvnaKk1KB3R1/+em8kgf4alRMRkXNzOqzs2LGDPn3MuYGdO3dWe0xbozdOpwvNs35OFkDncFh8Hfi6YOCj6snJg3v685c7I/Dx1n8jIiJyfk5/BH377beu6Ie4iSIr3PEpHDoNcUGwZLx5OKGzvvwxj2ffycQwYOSlATw8NRwvTwUVERG5MBUKSAWbHe7/AjanmgFlyXiIccEW+pv3FvGP/5hB5Yargrh/sk5OFhGR2nO6WGDevHm8/vrrZ11//fXXmT9/vrPNy0ViGPCX1fDVYfD1hH+Pc83pyclppTzx6klsdhgxIIAHblJQERGRunE6rLzyyit07tz5rOvdunXj5ZdfdrZ5uUgWboB3doAFWDAGLm3lfJvZeTYeeekkeYUG3dr58Kdbw1XHJCIideZ0WElLSyM2Nvas65GRkaSmpjrbvFwE/90Nz/5o3n5iKFzT0fk2S0oN/vLKKVJPWYkN9+Svv41UMa2IiDjE6bCSkJDAunXrzrq+bt064uLinG2+3v3vf/+jU6dOdOzYkddee62hu3PRrT4Cs742b/+uH9zWy/k2DcPg7//JYMehYgLLzvppEeyC0w5FRKRZcrrA9q677mLGjBmUlpYyfPhwAFatWsXDDz/MQw895HQH65PVamXmzJl8++23hISE0LdvXyZOnEhYWFhDd+2i2JEOv/scbAZM6AwPD3FNu+98kcPX6wvw8IAn7o6kday20BcREcc5HVYefvhhMjMzue+++ygpKQHAz8+PWbNm8cgjjzjdwfq0fv16unXrRqtWZoHGNddcw4oVK/jVr37VwD2rf/klMP1zKCiFyxPgmZGu2fTtm435vPG/bABm3BxGv85+zjcqIiLNmtPTQBaLhfnz53Py5El++ukntm3bRmZmJn/5y19c0b/zWrt2LePGjSMuLg6LxcLy5cvPes5LL71E27Zt8fPzo1+/fnz33XcVj6WkpFQEFYD4+HiOHz9e7/12B0+thV+yoVUwvHQt+LhglmbX4WLmv5UBwOQRwVx3uQvWPYuISLPnVFgpLS1l2LBh7N+/n6CgIAYMGED37t3x9XXBLmK1kJ+fT69evXjxxRdrfHzp0qXMmDGDxx57jC1btnDFFVcwduxYkpOTAbO24kznW61SXFxMTk5Ota/G6KtD8P4uc+XPc6Mh1AV/XKmnrMx++SSlVhjS0597JrRwvlERERGcDCve3t7s3LmzwZajjh07ljlz5jBx4sQaH3/uuee48847ueuuu+jSpQsLFiwgISGBRYsWAdCqVatqIynHjh2rcWVTuXnz5hEaGlrxlZCQ4Nof6CJIz4dZq8zbv+0Hl8U732ZeoZ1HF50kK89OhwRvHr09HE/tpSIiIi7i9DTQb37zG/7973+7oi8uVVJSwqZNmxg9enS166NHj+aHH34A4NJLL2Xnzp0cP36c3NxcPv/8c8aMGXPONh955BGys7Mrvo4ePVqvP4OrGQb8aSVkFkLXSJh5mfNtWm0GT756il9SSwkP9eTp30Xi76uDCUVExHWcLrAtKSnhtddeY+XKlfTv35/AwMBqjzfUqcunTp3CZrMRHR1d7Xp0dDRpaWkAeHl58Y9//INhw4Zht9t5+OGHCQ8/97atvr6+F22Kqz68vR1W/2LuUPvCGOcPJzQMg3/9v9Ns2luEn4+FufdFEtlCJziIiIhrOf3JsnPnTvr27QvA/v37qz3mDruVntkHwzCqXbv++uu5/vrrL3a3LrqDmfD09+btRy53zVb6H36by6ff5WGxwGN3hNMxwcf5RkVERM7QZE9djoiIwNPTs2IUpVx6evpZoy1NXYkNHlxhnqh8ZSJMc8HGbz9sL2DRh1kA3DuxBUN6BjjfqIiISA1cNma/e/dukpOTK/ZaAXNUY9y4ca56izrx8fGhX79+rFy5kgkTJlRcX7lyJTfccINTbS9cuJCFCxdis9mc7eZFseBn2JkOLfzg76Oc30/l4NES5ryRgWHAuMuDuHF4sGs6KiIiUgOnw8rhw4eZMGECO3bswGKxVCwHLp9qqc8P9Ly8PA4ePFhxPykpia1btxIWFkZiYiIzZ85k6tSp9O/fn0GDBrF48WKSk5O59957nXrf6dOnM336dHJycggNDXX2x6hXG47Doo3m7bnDIdrJrU9OZVl5dNFJiooN+nX24/4pLd1iuk9ERJoup5dtPPjgg7Rt25YTJ04QEBDArl27WLt2Lf3792f16tUu6OK5bdy4kT59+tCnTx8AZs6cSZ8+fSo2pJsyZQoLFizgqaeeonfv3qxdu5bPP/+c1q1b12u/3EVuMcz4CuwGTO4C1zp5QGGp1WD2y6c4lWWjdYwX/3dXBF6eCioiIlK/LEZNO6PVQUREBN988w09e/YkNDSU9evX06lTJ7755hseeughtmzZ4qq+up3ykZXs7GxCQkIaujtneegr+GAPxIfAl7dAsJMLmf79cRb/WZFDSKAHL82KIS5CK39ERMQxdfkMdXpkxWazERRkzi1ERESQkpICQOvWrdm3b5+zzYuDPjtgBhUPCywY7XxQ2Z1UzHtfmTv2zrwlTEFFREQuGqc/cbp378727dtp164dAwcO5JlnnsHHx4fFixfTrl07V/TR7bh7gW1aHjz6jXn7vv4woNX5n38hhcV2/rYkA7sBIy8N4Mo+WvkjIiIXj9PTQCtWrCA/P5+JEydy+PBhrrvuOvbu3Ut4eDhLly5l+PDhruqr23HHaSC7Ab9ZDt8lQ48oWHYTeDt5SOE/l2ayfE0eES08ef3xWIICtEOtiIg4py6foU6PrFTdnr5du3bs3r2bzMxMWrbUKpGG8PZ2M6j4ecGCMc4Hlc17i1i+Jg+AP90apqAiIiIXXb0UHoSFhdVHs3IBWUXwjx/N249cDh2c/GPIK7TzzNsZAFx/RRADuvo72UMREZG60/8mNyH//Bmyi6FzOEzt4Xx7C/97mvTTNuIivfjtxBbONygiIuIAhZUm4kgWvLXdvP3YFeDp5J/s91sLWPFTPhYL/Pk34TpJWUREGow+gRywcOFCunbtyoABAxq6KxX+tg5K7XBVa7jSyT3vTufaeO7dTACmjAyme/vGe9K0iIg0fgorDpg+fTq7d+9mw4YNDd0VwNxS/4uD5p4qj13uXFuGYfD8u5lk5dlpG+fNbde1cEkfRUREHKWw0sjZDZjznXl7SjfoFOFceyvXF/D9tkI8PeDP08Lx8daKLhERaVguCyuDBw8mOjraVc1JLX26H7aegEBvmHmZc22lZ1r51/8zp39+c20oHRN8XNBDERER57hs6fKUKVM4deqUq5qTWiiywjPrzNv39oeoQMfbMgyDZ9/JJL/QoHMbH24Z7R6b3ImIiLgsrDz44IOuakpq6Y2tcCwXYoLg7j7OtfXJ2jw27S3C19vCn6eF46nTlEVExE2oZqWRyiiAhWX1vQ8PBn9vx9s6ll7KK8uyALh7fAsSo51oTERExMUUVhzgDkuXn/8ZckugexRM6Ox4Oza7wd+WZFBUYtD7El/GXxXkuk6KiIi4gMKKAxp66fKBTHh3h3n7scvNJcuOWroyh91JJQT4WXh4ajgezjQmIiJSD+o1rOTk5NRn883WvO/BZsDItjA4wfF2Dh0r4c3/ZQMwfXJLYsLr5agoERERpzgcVv7+97+f9/GcnBxGjx7taPNyDuuOwqok8LSYhxU6yjAMFryfidUGg3r4c/VlTiwlEhERqUcOh5XZs2fzxhtv1PhYXl4eY8aM0ciKi9nslRvA/bqHc6cqr9lSyK7DJfh6W5hxc0ssFk3/iIiIe3I4rLz99tvcd999LF++vNr1vLw8Ro8eTWZmJt9++62z/ZMqlu2F3Sch2AdmDHS8nZJSg1eXnQZgyqhgIltq+kdERNyXw59SN954I1lZWdxyyy189tlnDBs2jLy8PK6++mpOnTrFmjVrtKOtCxWWwrM/mrd/PwDCAxxva9nqXFIzbISHejJllDZ/ExER9+bU/1LfddddZGZmMn78eD7++GNmz55NWloaa9asITY21lV9FOD/7Ya0PIgPhtt6O95OVq6Nd74wi2rvvD4Uf18tCBMREffm9Pj/ww8/zOnTpxkxYgRt2rRhzZo1tGrVyhV9kzI2O7y2xbx9Tz/wc+JPbcnn2eQXGXRI8Gb0QBXVioiI+3P4Y2/ixInV7nt7exMREcEDDzxQ7fpHH33k6Fu4rYULF7Jw4UJsNttFeb8VhyA5G1r4weSujrfzS2opn36XB8DvJrbUnioiItIoOBxWQkNDq93/1a9+5XRnGovp06czffp0cnJyzvo9uJphwOLN5u1be0CAEzvhv7LsNHY7DO7pT59Ofq7poIiISD1zOKyca9myuNbGVNiSBj6eMK2XE+3sKeSnnUV4esA9E1q4rH8iIiL1zenqyqNHj7qiH3IOr5aNqkzoDFEOlpjY7AYvf5gFwA1XBumgQhERaVScLrBt3bo1LVu2pFevXvTq1YvevXvTq1cviouLWbhwIW+99ZYr+tksJZ2Grw6Zt+/q43g7X/6Yz+GUUoL8LUy9pn6nrURERFzN6bBy+PBhtm7dytatW9myZQsffPABKSkpAISEaA8PZ7y2BQxgeBu4JNyxNgqK7LzxaRYAU68JJTTI01XdExERuSicDitt2rShTZs2jB8/vuLajz/+yLRp05g/f76zzTdbGQXw393m7Xv6Od7O+ytzyMyx0yrSi/FXBbumcyIiIhdRvewINmjQIF544QXmzJlTH803C+/sgGIb9IiCyxzctiY908r/+zoXMItqvb20VFlERBofp8NKaWlpjdc7duzIrl27nG2+WSqywpJt5u17+oKjZwy+9kkWJaUGPTv4cnkvf9d1UERE5CJyehooMDCQrl270qdPH3r37k2fPn2Ii4vjX//6F6NHj3ZFH5udj/ZARiG0CoZrOjrWxt4jxXy9vgCA+27UqcoiItJ4OR1WvvnmG7Zt28a2bdv4z3/+w6OPPkphYSEAo0eP5rHHHqNnz5707NmTLl26ON1hd1CfO9jajcqt9e/oDV4OjH0ZhsGisqXKowcGckmij8v6JyIicrFZDMMwXNmg3W5n3759FSuEyoNMenr6Rdue/mIp38E2OzvbZSufvj4Md34KIT7w450Q5EDOWLulgCdePYWvt4W3noglsqXTmVRERMSl6vIZ6vJPMQ8PD7p06UKXLl2qbcF/4sQJV79Vk1S+tf4tPRwLKiWlBq8sywLgppHBCioiItLoOVRgu337dux2e62fv2vXLsLDHdwopBnZmgY/Hzenfm7v7VgbH6/NJfWUlbAQD24epX1uRESk8XMorPTp04eMjIxaP3/QoEEkJyc78lbNSvnW+jd0gpigur++pNSoWKp823Ut8Perl5XpIiIiF5VDcwSGYTB79mwCAgJq9fySkhJH3qZZSc6Gzw+at+92cGv9levzyci2EdHCkzGXOXiQkIiIiJtxKKxceeWV7Nu3r9bPHzRoEP7+2ufjfN7cZq4EujIRukTW/fU2u8HSlTkATB4RrA3gRESkyXAorKxevdrF3WjeSmywbK95+7ZejrWxblshx9KtBPlbuHaIA3NIIiIibkpFDW5g9RHILITIALiqTd1fbxgG731ljqrccFUwAapVERGRJkSfam7ggz3m9/GdHdsEbuv+Yvb9UoKPt4WJw3RYoYiINC0KKw3sdCF8k2TevtHBDX7fL6tVGTsokJbBni7qmYiIiHtQWGlgH++HUjt0i4TOEXV//cGjJWzYXYSHBSaP1L4qIiLS9Lg0rBQUFJzzsaSkJFe+VYNauHAhXbt2ZcCAAU639eFu87ujoyrvlY2qDO0XQFyEdqsVEZGmx6VhJTQ0lMcffxyr1XrWY5MmTXLlWzWo6dOns3v3bjZs2OBUO/szYHu6WadyQ6e6vz7llJU1m8yAqN1qRUSkqXJpWGnfvj0pKSkMHDiQvXv3VnvMxeclNgnlhbXD20B47fbXq+a/X+dgN2BAVz86JOhkZRERaZpcGlYCAgJ4/fXXmT17NqNHj+aFF16oeMxi0SZlVVntlXurTHJgCuh0ro0vfswH4FejNaoiIiJNV70U2I4fP54NGzawatUqRo4cyfHjx+vjbRq175MhPR9a+sHwtnV//Uff5lJSatC5jQ+9Ovq6voMiIiJuwqUVmVWneqKjo/nkk09YvHgxl112GYWFha58q0bvw7IpoOs7gU8dVxsXFNn5eI15YOHNo0I0aiUiIk2aS0dW7r777rOu3XPPPaxZs4Ybb7zRlW/VqOUUw4pD5m1HVgF9ti6PvEKDhGgvLu+lM5dERKRpc+nIyrXXXkt2djahoaEAfP3113zyySe0bt26Wv1Kc/fZASi2Qccw6BFVt9eWWg3+u8ocVZkyMgQPD42qiIhI0+bSkZXJkyeTn28WfW7cuJGbbrqJxMREdu7cyb333uvKt2rU/lu+t0pXqOsMztcb8jmVZSM81JORlwa6vnMiIiJuxqUjK0VFRcTFxQHwzjvvcM899/DHP/4RwzDo2bOnK9+q0Uo6DZtSwcMCEzrX7bV2u8HSsgMLJw0PxsdboyoiItL0uXRkxTCMiiLblStXMmrUKEDLlqv6sGy58hWJEF3HgZEfdxSSfMJKoL+FcZcHub5zIiIibsilIys33XQTo0aNomXLlnh6ejJs2DAADh8+THCwTgO2G/BR2SqguhbWGobBe2WjKjdcGUygv451EhGR5sGlYWX27NmMHDmStLQ0Ro0ahYeH+YFqtVp58cUXXflWjdJPx+B4LoT4wOj2dXvt7qQSdieV4O0FE4cp+ImISPPh8v89HzRoEH369CEoqHKa4pJLLqnxvKDmpnxvlesuAb86xsT/fZ8HwPD+gYSF1HFjFhERkUasXuYSJk+ezC+//FJx/6uvvuKee+6pj7dqNIqs8PlB83Zdt9fPK7SzZrN5YOG1Q1SrIiIizUu9hJXXX3+dG2+8kSNHjrB06VIef/xxVqxYUR9v1WisOwoFpRAbBP1i6/babzbkU1Ri0DrWm27tdGChiIg0Ly6tWSnXo0cPXn/9da699lqioqL4+uuvCQlp3oftrTxsfh/Zru57q3y2zpwCumZwoFZWiYhIs+NwWFmxYkW1IlqAAQMGVPswzcrKwmKxMHLkSADWr1/vRFcbL7sBX5eFlTHt6vba/cklHDhaircXjB6oTeBERKT5cTisXHPNNaSmphIVVblf/AcffOCSTjU1W9LgZAEE+8DA+Lq99vOyUZXLewcQGqTCWhERaX4cDitVT1gu17p1a6c601gsXLiQhQsXYrPZavX8lWWHFg5rU7cTlguL7azaYB5fcO1gFdaKiEjzpJ3FHDB9+nR2797Nhg0bavX8r8qmgEbVcQpo7ZYC8osMYiO86H2Jbx17KSIi0jQ4FVZefPFFVqxYwalTp1zVnybn0Gnzy9sDhrap22s/W2eOqlwzOFCnK4uISLPl1GqgRYsWMWfOHCwWC61ataJv377069ePvn370rdvX2Jj67hGtwkqnwIaFA8hdRgcOZJays5DxXh4wNWDNAUkIiLNl1NhZdeuXVitVrZs2cLmzZvZvHkzr732GkePHsVisRAdHU1KSoqr+tooOToF9MUPZmHtoO7+hIeqsFZERJovh8NK+RLluLg44uLiuPbaaysey8zMZOPGjWzdutXpDjZmJ/Nhc6p5uy5hpaTUYMVPZVNA2rFWRESaOZeuBioXFhbG6NGjGT16tKPNNwlfJ4EB9IyC2DqcPbhuWwE5+XYiWnhyaVe/euufiIhIY+Bwge0XX3xBaGioK/vS5JTvWjuqjicsf/6DOaoydlAgnp4qrBURkebN4ZGVMWPGuLIfTU5BKXyfbN4eXYcpoJRTVjbtLcJigbHaW0VERET7rNSXtb9AsQ0SQqBTeO1f90XZjrX9OvsRE14vRzeJiIg0Kgor9aR8FdDo9rU/uNBmM/hShbUiIiLVKKzUA6sdViWZt+syBfTTrkIysm20CPJgSE//+umciIhII6OwUg82pkBWEbTwg/5xtX/dZ9+bU0CjLwvE20uFtSIiIqCwUi/Kp4BGtgWvWv6GT2ZZWb+rCIBrVFgrIiJSQWHFxQwDvirbYr8uG8Gt+DEfuwE9OviSGONdP50TERFphBRWXGxfBhzNAV9PuLJ17V/3zcYCAK4eFFhPPRMREWmcFFZcbPUR8/uQBAio5QDJkdRSjqSW4uUJV/QKqLe+iYiINEYKKy627qj5/YrE2r9mzWZzVKVfZz+CAvRHIiIiUpU+GV2o2Arryw6ZHpJQ+9eVh5Wh/TSqIiIiciaFFRfakgZFVogMgEtquWvtL1WmgIb0VFgRERE5k8KKC5VPAQ2Kr/2utWu2aApIRETkfPTp6ELlYcWRKaCr+mpURUREpCYKKy6SVwLbTpi3h9SyuPaX1FKSUsqmgLQKSEREpEYKKy7y83HzTKDEUPOk5doonwLq29mPYE0BiYiI1EifkC7izBTQUE0BiYiInFOzDysTJkygZcuW3HjjjU6180Mdw0pymqaAREREaqPZh5UHHniAt956y6k2Mgpgzynz9qD42r2mfFRFU0AiIiLn1+w/JYcNG0ZwcLBTbZRvBNc5HCJqOUiiVUAiIiK149ZhZe3atYwbN464uDgsFgvLly8/6zkvvfQSbdu2xc/Pj379+vHdd99d9H7+dMz8XpcpoMMppXh6wJCe/vXXMRERkSbArcNKfn4+vXr14sUXX6zx8aVLlzJjxgwee+wxtmzZwhVXXMHYsWNJTk6ueE6/fv3o3r37WV8pKSku6+dPx83vtQ0rVc8CCgn0dFk/REREmiKvhu7A+YwdO5axY8ee8/HnnnuOO++8k7vuuguABQsWsGLFChYtWsS8efMA2LRpk8v6U1xcTHFxccX9nJwcAI7lgLc/XNqqdu1UTAHpLCAREZELcuuRlfMpKSlh06ZNjB49utr10aNH88MPP9TLe86bN4/Q0NCKr4SEyqGUXjEQ7HvhNpJPaApIRESkLhptWDl16hQ2m43o6Ohq16Ojo0lLS6t1O2PGjGHy5Ml8/vnnxMfHs2HDhnM+95FHHiE7O7vi6+jRoxWPaQpIRESkfrj1NFBtWM44MdAwjLOunc+KFStq/VxfX198fWsePqlrWNEqIBERkdpptCMrEREReHp6njWKkp6eftZoS33z9YS+MRd+3tETpRw+XjYF1EtTQCIiIrXRaMOKj48P/fr1Y+XKldWur1y5ksGDB1/UvvSLBd9ajFFV3QhOU0AiIiK149bTQHl5eRw8eLDiflJSElu3biUsLIzExERmzpzJ1KlT6d+/P4MGDWLx4sUkJydz77331mu/Fi5cyMKFC7HZbAAMrO0qoC06C0hERKSuLIZhGA3diXNZvXo1w4YNO+v6tGnTePPNNwFzU7hnnnmG1NRUunfvzvPPP8+VV155UfqXk5NDaGgo6w5kM7jD+Y9aTs+0cvPjKXhY4MP5rQgN0siKiIg0X+WfodnZ2YSEnP8z1K3Dirsr/0Vnns6mZYvz/6L/930ez72bSbd2Pvzrj7UocBEREWnC6hJWGm3NijvxrMVv8eddhQBc2lWFtSIiInWhsHIRlFoNtuwrAuDSbn4N3BsREZHGRWHFAQsXLqRr164MGDCgVs/fdbiYgiKDFkEedEzwqefeiYiINC0KKw6YPn06u3fvPu9ut1WtL5sC6t/VDw+P2m9YJyIiIgorF8X6XeYU0MBuqlcRERGpK4WVenYyy8rhlFIsFujfRfUqIiIidaWwUs82lI2qdG7to71VREREHKCwUs/W7zbrVQZ01aiKiIiIIxRWHFDb1UBWm8GmPapXERERcYbCigNquxpod1Ix+UUGIYEeXNJaS5ZFREQcobBSj8pXAQ3o6oenliyLiIg4RGGlHpXvrzJAW+yLiIg4TGGlnmRk2zh4rBRQca2IiIgzFFbqyYayVUCdEn1oGawlyyIiIo5SWKknFfUqOrhQRETEKQorDrjQ0mWbzWDT3rJTllWvIiIi4hSFFQdcaOnyniMl5BbYCQ7woEtbLVkWERFxhsJKPSjftbZ/Fy1ZFhERcZbCSj2our+KiIiIOEdhxcVO59rYn1wCqF5FRETEFRRWXKz8LKAO8d6EhWrJsoiIiLMUVlxs6wEzrPTtrCkgERERV1BYcbHtB4oB6NVRYUVERMQVFFYccK59VjKybRxLt2KxQI/2vg3UOxERkaZFYcUB59pnZVvZFFD7eG+CAvSrFRERcQV9orrQNk0BiYiIuJzCigttLxtZ6dVBU0AiIiKuorDiIqdzbfySZgWgh8KKiIiIyyisuEj5KqB2cd6EBml/FREREVdRWHGR8uLanh01qiIiIuJKCisuov1VRERE6ofCigtk59k4nFIKaGRFRETE1RRWXGD3YXNUpXWMFy2DVa8iIiLiSgorDjhzB9sdZWGlp6aAREREXE5hxQFn7mC781B5vYqmgERERFxNYcUFksrqVVRcKyIi4noKKy5gGBAf5UV4qOpVREREXE1hxUW0CkhERKR+KKy4SK8OmgISERGpDworLqLiWhERkfqhsOICMWGeRIV5NXQ3REREmiSFFRfo1l5TQCIiIvVFYcUFurf3aeguiIiINFkKKy7QQyMrIiIi9UZhxQWiw7S/ioiISH1RWHEBi8XS0F0QERFpshRWHHDmQYYiIiJSfyyGYRgN3YnGKicnh9DQULKzswkJCWno7oiIiDQadfkM1ciKiIiIuDWFFREREXFrCisiIiLi1hRWRERExK0prIiIiIhbU1gRERERt6awIiIiIm5NYUVERETcmldDd6AxK99PLycnp4F7IiIi0riUf3bWZm9ahRUn5ObmApCQkNDAPREREWmccnNzCQ0NPe9ztN2+E+x2OykpKQQHB7vdYYY5OTkkJCRw9OhRHQVQS/qdOUa/t7rT78wx+r3VnTv/zgzDIDc3l7i4ODw8zl+VopEVJ3h4eBAfH9/Q3TivkJAQt/sP1N3pd+YY/d7qTr8zx+j3Vnfu+ju70IhKORXYioiIiFtTWBERERG3prDSRPn6+vJ///d/+Pr6NnRXGg39zhyj31vd6XfmGP3e6q6p/M5UYCsiIiJuTSMrIiIi4tYUVkRERMStKayIiIiIW1NYEREREbemsCIiIiJuTWGliVm7di3jxo0jLi4Oi8XC8uXLG7pLbm/evHkMGDCA4OBgoqKiGD9+PPv27Wvobrm1RYsW0bNnz4pdMQcNGsQXX3zR0N1qdObNm4fFYmHGjBkN3RW39cQTT2CxWKp9xcTENHS3GoXjx49z6623Eh4eTkBAAL1792bTpk0N3S2HKKw0Mfn5+fTq1YsXX3yxobvSaKxZs4bp06fz008/sXLlSqxWK6NHjyY/P7+hu+a24uPj+dvf/sbGjRvZuHEjw4cP54YbbmDXrl0N3bVGY8OGDSxevJiePXs2dFfcXrdu3UhNTa342rFjR0N3ye2dPn2aIUOG4O3tzRdffMHu3bv5xz/+QYsWLRq6aw7R2UBNzNixYxk7dmxDd6NR+fLLL6vdf+ONN4iKimLTpk1ceeWVDdQr9zZu3Lhq959++mkWLVrETz/9RLdu3RqoV41HXl4ev/71r3n11VeZM2dOQ3fH7Xl5eWk0pY7mz59PQkICb7zxRsW1Nm3aNFyHnKSRFZEzZGdnAxAWFtbAPWkcbDYb77//Pvn5+QwaNKihu9MoTJ8+nWuvvZaRI0c2dFcahQMHDhAXF0fbtm25+eabOXz4cEN3ye198skn9O/fn8mTJxMVFUWfPn149dVXG7pbDlNYEanCMAxmzpzJ5ZdfTvfu3Ru6O25tx44dBAUF4evry7333suyZcvo2rVrQ3fL7b3//vts3ryZefPmNXRXGoWBAwfy1ltvsWLFCl599VXS0tIYPHgwGRkZDd01t3b48GEWLVpEx44dWbFiBffeey8PPPAAb731VkN3zSGaBhKp4ve//z3bt2/n+++/b+iuuL1OnTqxdetWsrKy+PDDD5k2bRpr1qxRYDmPo0eP8uCDD/LVV1/h5+fX0N1pFKpOa/fo0YNBgwbRvn17lixZwsyZMxuwZ+7NbrfTv39/5s6dC0CfPn3YtWsXixYt4je/+U0D967uNLIiUub+++/nk08+4dtvvyU+Pr6hu+P2fHx86NChA/3792fevHn06tWLF154oaG75dY2bdpEeno6/fr1w8vLCy8vL9asWcM///lPvLy8sNlsDd1FtxcYGEiPHj04cOBAQ3fFrcXGxp71Pw5dunQhOTm5gXrkHI2sSLNnGAb3338/y5YtY/Xq1bRt27ahu9QoGYZBcXFxQ3fDrY0YMeKslSy33347nTt3ZtasWXh6ejZQzxqP4uJi9uzZwxVXXNHQXXFrQ4YMOWsLhv3799O6desG6pFzFFaamLy8PA4ePFhxPykpia1btxIWFkZiYmID9sx9TZ8+nXfffZePP/6Y4OBg0tLSAAgNDcXf37+Be+eeHn30UcaOHUtCQgK5ubm8//77rF69+qyVVVJdcHDwWbVQgYGBhIeHq0bqHP74xz8ybtw4EhMTSU9PZ86cOeTk5DBt2rSG7ppb+8Mf/sDgwYOZO3cuN910E+vXr2fx4sUsXry4obvmGEOalG+//dYAzvqaNm1aQ3fNbdX0+wKMN954o6G75rbuuOMOo3Xr1oaPj48RGRlpjBgxwvjqq68auluN0lVXXWU8+OCDDd0NtzVlyhQjNjbW8Pb2NuLi4oyJEycau3btauhuNQqffvqp0b17d8PX19fo3LmzsXjx4obuksMshmEYDZSTRERERC5IBbYiIiLi1hRWRERExK0prIiIiIhbU1gRERERt6awIiIiIm5NYUVERETcmsKKiIiIuDWFFREREXFrCisiIiLi1hRWRKTevPDCC7Rt25aAgADGjx9PdnZ2jc8bOnQoFosFi8XC1q1bL24na2no0KHMmDGjTq+57bbbKn6u5cuX10u/RJoDhRURqRePPvooL774IkuWLOH7779ny5YtPPnkk+d8/t13301qaqrbHuj30Ucf8de//rVOr3nhhRdITU2tpx6JNB8KKyLichs2bGD+/PksXbqUK6+8kr59+/Lb3/6W//3vf+d8TUBAADExMXh5OX4YfElJicOvvZCwsDCCg4Pr9JrQ0FBiYmLqqUcizYfCioi43N///neGDx9O3759K65FRkZy6tSpOrXz5Zdfcvnll9OiRQvCw8O57rrrOHToUMXjQ4cO5fe//z0zZ84kIiKCUaNGAWC325k/fz4dOnTA19eXxMREnn766Wqvu//++5kxYwYtW7YkOjqaxYsXk5+fz+23305wcDDt27fniy++qPaaqtNAQ4cO5YEHHuDhhx8mLCyMmJgYnnjiiTr+pkSkNhRWRMSliouL+fTTT5kwYUK164WFhYSGhtaprfz8fGbOnMmGDRtYtWoVHh4eTJgwAbvdXvGcJUuW4OXlxbp163jllVcAeOSRR5g/fz6zZ89m9+7dvPvuu0RHR1dre8mSJURERLB+/Xruv/9+fve73zF58mQGDx7M5s2bGTNmDFOnTqWgoOCc/VuyZAmBgYH8/PPPPPPMMzz11FOsXLmyTj+jiNSCISLiQj/88IMBGH5+fkZgYGDFl4+PjzFmzJgaX3PVVVcZDz744AXbTk9PNwBjx44dFa/r3bt3tefk5OQYvr6+xquvvnrOdq666irj8ssvr7hvtVqNwMBAY+rUqRXXUlNTDcD48ccfa+zjmW0YhmEMGDDAmDVr1lnvBxjLli274M8nIjVzfHJYRKQG+/fvx8/Pjx07dlS7fv311zNkyJA6tXXo0CFmz57NTz/9xKlTpypGVJKTkysKcfv371/tNXv27KG4uJgRI0act+2ePXtW3Pb09CQ8PJwePXpUXCsfiUlPT69VGwCxsbHnfb6IOEZhRURcKicnh6ioKDp06FBxLTk5mb179zJp0qQ6tTVu3DgSEhJ49dVXiYuLw263071792qFtIGBgdVe4+/vX6u2vb29q923WCzVrlksFoBqU061aeN8zxcRx6hmRURcKiIigpycHAzDqLj29NNPc80119C1a9dat5ORkcGePXt4/PHHGTFiBF26dOH06dMXfF3Hjh3x9/dn1apVDvVfRNyPRlZExKWGDx9OUVERf/vb3/jVr37Fu+++yyeffML69evr1E7Lli0JDw9n8eLFxMbGkpyczJ///OcLvs7Pz49Zs2bx8MMP4+Pjw5AhQzh58iS7du3izjvvdPTHEpEGpJEVEXGp6Oho3nzzTRYtWkTXrl354Ycf+P7770lISKhTOx4eHrz//vts2rSJ7t2784c//IFnn322Vq+dPXs2Dz30EH/5y1/o0qULU6ZMUS2JSCNmMaqO1YqINIChQ4fSu3dvFixY0NBdqRcWi4Vly5Yxfvz4hu6KSKOkkRURcQsvvfQSQUFBZ60iaszuvfdegoKCGrobIo2eRlZEpMEdP36cwsJCABITE/Hx8WngHrlGeno6OTk5gLms+cyVSyJSOworIiIi4tY0DSQiIiJuTWFFRERE3JrCioiIiLg1hRURERFxaworIiIi4tYUVkRERMStKayIiIiIW1NYEREREbemsCIiIiJu7f8DbxHK2O5tiEUAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "f, ax = plt.subplots(1, 1, figsize=(6, 4))\n", "\n", "ax.semilogy(\n", " theta, T_kSZ_150 * glasz.constants.sr2sqarcmin, label=\"150 GHz\", color=\"dodgerblue\"\n", ")\n", "ax.semilogy(\n", " theta, T_kSZ_090 * glasz.constants.sr2sqarcmin, label=\"90 GHz\", color=\"royalblue\"\n", ")\n", "\n", "ax.set_xlabel(r\"$\\theta$ [arcmin]\")\n", "ax.set_ylabel(r\"$T_{\\rm kSZ}$ [$\\mu$K $\\cdot$ arcmin$^2$]\")\n", "ax.legend()\n", "ax.set_xlim(0.5, 6.5)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "glasz", "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.11" } }, "nbformat": 4, "nbformat_minor": 2 }