{ "cells": [ { "cell_type": "markdown", "id": "60cd9e50", "metadata": {}, "source": [ "# On the importance of having the correct survey dimension and k-sampling for intensity mapping\n", "\n", "`meer21cm` naturally deals with survey geometry, and therefore it is easier than most tools to have realistic forecasts of your desired signal with the correct\n", "lightcone dimensions and sampling.\n", "\n", "In particular, the clustering signal is anisotropic, due to signal properties such as RSD as well as observational effects such as the beam. Espeically for intensity mapping, where the transverse and line-of-sight direction have significantly different scales and resolutions, having the correct box lengths, grids, and therefore the correct k-sampling is crucial for forecasts.\n", "\n", "If you want to forecast the measurements of some power spectrum statistics, we highly recommend that you \n", "\n", "- **use `meer21cm` to calculate the survey lightcone and k-modes**\n", "- **build your model at the level of 3D k-space using your own tool** (since `meer21cm` only provides some simple model routine)\n", "- **apply window convolution, observation effects either yourself or using `meer21cm`**\n", "- **and finally average into the (usually 1D) statistics you want to use for forecasts**.\n", "\n", "For this notebook, we will demonstrate the effect of simply having a non-cubic box and how important it is to have the correct k-modes at the 3D level.\n", "\n" ] }, { "cell_type": "markdown", "id": "36c04e1f", "metadata": {}, "source": [ "## Survey set-up\n", "\n", "Suppose you want to forecast for a MeerKLASS UHF-like survey, with ~7500 sqdeg for the southern patch from z-0.5 - 1.5:" ] }, { "cell_type": "code", "execution_count": 1, "id": "b8787629", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/zhaotingchen/miniconda3/envs/meer21cm/lib/python3.10/site-packages/halomod/halo_model.py:32: UserWarning: Warning: Some Halo-Exclusion models have significant speedup when using Numba\n", " from .halo_exclusion import Exclusion, NoExclusion\n" ] } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from meer21cm import util\n", "from meer21cm.plot import plot_map\n", "from meer21cm import PowerSpectrum" ] }, { "cell_type": "code", "execution_count": 2, "id": "0b6bf99c", "metadata": {}, "outputs": [], "source": [ "z_min = 0.5\n", "z_max = 1.5\n", "# frequency resolution in Hz, use your own number\n", "freq_resol = 0.5e6\n", "num_ch = int((util.redshift_to_freq(z_min)-util.redshift_to_freq(z_max))/freq_resol)\n", "# obtain the frequency channels\n", "nu = np.linspace(0,num_ch-1,num_ch) * freq_resol + util.redshift_to_freq(z_max)" ] }, { "cell_type": "code", "execution_count": 3, "id": "3f7c0102", "metadata": {}, "outputs": [], "source": [ "ra_range = [-45, 80]\n", "dec_range = [-70, 5]\n", "ang_resol = 1.0 # in deg, similar to beam size here for simplicity\n", "num_pix_x = 102 # one extra on each side as buffer\n", "num_pix_y = 102 # one extra on each side as buffer\n", "wcs, num_pix_x, num_pix_y = util.create_wcs_with_range(\n", " ra_range,\n", " dec_range,\n", " [ang_resol,ang_resol],\n", ")\n", "\n", "ps = PowerSpectrum(\n", " wproj=wcs,\n", " num_pix_x=num_pix_x,\n", " num_pix_y=num_pix_y,\n", " nu=nu,\n", " ra_range=ra_range,\n", " dec_range=dec_range,\n", ")" ] }, { "cell_type": "code", "execution_count": 4, "id": "b0d0c8d6", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAFqCAYAAAD87JuUAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAbsFJREFUeJztnQd0XPd15j+i9947SBAkCrvYe++dIinJG1mxvV7H9tnEceIWW1mv2zou8SYuWa8lO16JYu+9917BAhKF6L23ATAAZs+9EGB2DoCZee3+jp6GmBnM/AG8ee97t3x3kMVisUAQBEEQBEGnOCm9AEEQBEEQBHsiYkcQBEEQBF0jYkcQBEEQBF0jYkcQBEEQBF0jYkcQBEEQBF0jYkcQBEEQBF0jYkcQBEEQBF0jYkcQBEEQBF3jAoPT1dWFkpIS+Pr6YtCgQUovRxAEQRAEKyBP5MbGRkRFRcHJ6dWxG8OLHRI6sbGxSi9DEARBEIR+UFhYiJiYGHWnsX77299i5MiR8PPz423y5Mk4ePBg7+Otra348pe/jODgYPj4+GDt2rUoLy9/6jX27NmD5ORkDBs2DPv27evT+1NEp5tBn/46nt+cnFwxatRobNy4EW+8MQGuru4vfa5s+tw8Pb2xZMlSxdchm2zWbuvWvan4GmQzxrZ+/Xq7vG5ERBSWL1/BW3h45AueM+iZ87iKIzukxn7yk59g6NChHJL605/+hJUrV+LmzZtIS0vD3/3d32H//v3YunUr/P398ZWvfAVr1qzB+fPn+fvb2tpYDH344Yf8/X/913+NBQsWwM3Nzar3/0vqim6fTmO5u7th3LhxSEoajDt3MrB163Z0dnY98XzBKLi6un26r8jfXdAGrq6usr8KDqG5uQX+/gGor6+36euWlZVj7979CAgIwOTJEzkgcvHiJRQUFD7xLItVJSiKi53ly5c/9fUPf/hDjvZcunSJhdAf/vAHfPzxx5gzZw4/TqImJSWFH580aRKLHWdnZ4wePZofd3Fx4fusFTsvwtvbGxMnTkBkZASuX7+BP//5ImRcqrFxdXWB2WxWehmCIAiqo7S0FFFRkTYXOz3U1dXh4MHDfG6ePHkSpk2bgsuXryInJ8fq11Bc7DxJZ2cnR3Cam5s5nXX9+nU+wcybN6/3OcOHD0dcXBwuXrzIYoeU3nvvvYfIyEhWdz/4wQ+sCmm9CEqTTZ06GYGBgaweT5w4acOfTtAyJKI7OjqVXoYgCILqKCkpQ2pqCh48yLTr+5A2OHbsODw9PTBhwgSMGzcGW7Zs0Y7YycjIYHFD9TkkOHbu3InU1FTcunWLIzQUwnqS8PBwlJWV9X79/vvv42//9m+5Gru/QmfevDkIDg7B+fMXUFRUPOCfSdCj2JHIjiAIwrNUVlYgLGwGHIXJ1IrTp898mqqFdsQOFRaTsKEQ2LZt2/Duu+/i9OnTfXoNqufpK0ePHsUXv/hF/vf9+5koKSnt82sIxoA+VGZzh9LLEARBUB2dnV2vbf22B2Zzu9XPdfzqXgBFb5KSkrgY+Mc//jFGjRqFX/3qV4iIiEB7ezvn656EurHosf5y+PBhDBkyhAuZewqbqAVdEF6GRHYEQRBeTkNDA5eVqBVViJ0XGf1RkTGJH7qiPn78eO9jDx8+REFBAae9+gOlypYtW8Ynr5MnT3LXlyC8DonsCIIgvJzKykqEhYVCrSiexvrWt76FxYsXc9ExOSFS59WpU6c4+kKpqc997nP42te+hqCgIFaNX/3qV1noUHFyf/Dw8ODIEL1ejxoVBOsiOyJ2BEEQXkRFRSVnXLKzre+QMpTYqaiowF/91V9x6xqJGzIYJKEzf/58fvyXv/wl5wLJTJCiPQsXLsRvfvObAb1nj9ARBGshewPa/wRBEIQXR3ZGjhwBtaK42CEfnddFYn7961/zJghK4eQ0iNOrgiAIwvM0NTWzD45aUWXNjiCoDSpkJ4duQRAEQXuI2BEEKxg0yEnEjiAIwisg+5j+2MA4AhE7gmBlZEfSWIIgCC+nsrIKISHBUCMidgTBypodiewIgiC8nNraWgQFBUKNiNgRBCuQmh1BEIRXU1NTy7Ml1YiIHUGwAqnZEQRBeDU07eDZWZZqQcSOIFiB1OwIgiC8ms7OTri4OEONiNgRBCvomaEmCIIgaA8RO4JgBRZLlwgeQRCE19DY2ARfXx+oDRE7gmAFXV0WHlsiCIIgaK9uR47egmAFVK8jkR1BEIRX09DQCF9fX6gNETuCYGUaSyI7giAIr6ahoQF+fn5QG3L0FgQrIzsidgRBEF5NY2Mj/PwksiMImkRqdgRBEF6PpLEEQcNIzY4gCMLraW9vh5ubG9SGiB1BsAJyT5bIjiAIgjaRo7cgWIHU7AiCIGgXF6UXIAhaoKOjA56eHjAS5o5jSi9BeAmuLvOUXoIgvBSz2QxXVxeYzR1QC3KpKghWfnhdXFyVXoYgCILqMZlM8PT0gpoQsSMIVl+piNgRBEF4HS0tJnh5eUJNSBpLEPoQllUDkl4SXrcPUEH91i338J9//sKA3kfSZUJ/aGlpgaenusSORHYEwcqaHYnsCFqho6MLLi5yeBeUS2N5eUkaSxA0h6SxBK2JHVdXObwLykV2vCSNJQjag7oKXFz693GRtJPgaMzmLri4Og/8dWy870pazDiRnaioKKgJkf6CYAUS2RG0hNncCVdJYwkK0damPhdl+TQIghV0dIjYETRWsyNpLEEhzGYSO+o6XsqnQRCsQAaBCpqL7NggjSUI/aG93ay6yI7U7AhCH5D6G0ErNTtqTGP15/MjdT7aHAbqqrJIuPo+DYIgCMKAkDSWoCSdnZ39buiwF/JpEARB0BmSxhKEp1GX9BIEFYfUt2y+a7e1CIKt01jubs6GTh1L+kt4EonsCIKVODkNQmdnl9LLEKyERibQ34s2KjCnr2kzApLGEoSnkciOIFiJu4cL2to64eUlJ5H+QEKDTsKtrR1oa+tAW2snWum2rePT+zrR3t6Jjo5Ojkx0mLs4HUPf02WxYBAG/eW18PTXL2QQ4DRoEN+yyKH/XqF1el6Tbru//fnXp5ejMQyUInJzd4abmzNHUNzcXT697b7Pzc0F7u7OcHd3gaenC5ydHbvPSBpLEJ5GxI6geRzVIeXh7swnZS8vdXUZKAGJBxImLS3mF2ztaGk2w9Ta8VQkhcQDRRs8SBh4dAsB2jw8XODj44aQYBe4ujnzmANXF2d+Lp2wSVxQVE0NUISIxBeJifa2TrS1k0Dr4N8FibXmJjNq202f3t+JttZuIdfZ+bTKop/R09OVhZCnF926wtvblX8PPj7ufP8gUlYDcVBWYTeWI5HOL+FJROwIQh8iO3Ti0juU9mlqbEdjYxsamz69pa2hnaMwPQEPN1dnFn68ebvxbUCAR+99dAJXi0ixFfTzdEdunOHt3b/XIAFIYsRkMndvLR18W13Vgvy8OjQ1tcNk6uiNMFF0ytvbDd4shNzg6+sGf38P+Pm5s1B8kSgiMUZrFAShGxE7gmAlFJHgk73GochEQ0Mr6mpbUVfXivr67tt2cyc/7uzk1HtS9fV1R2iINxITA/nflJoZSMRBoFTYXwQTiRZroknNze0sgmhrbGhDSXEj/90oDdjzmr4+n4ogf3dUV7fA3N4dWZO/lyCI2BEEXUZ2KIVCJ7yaGhPf1taYek+MVD/i7++OgABPvo2J9eeIjEQC1BtNIqFJ26sEUVNjG+ob2lgE1de14vqNUpjOFbLgcXZxQmCABwIDPREY5ImgIE+ODIkQEoyCiB1B1ajJsZhSMw0NbVALdIKjExuLmupuUdPQ2MaFuFT7Qie04GAvDBkShKDxnpxWEvQriPw4quOB2Fh/3L9fiaVLh/WmESmaR9G7mpoWVFQ0IzOziiNELIScnXhfCQnxQkioN98aVfhae7yR2h7tIWJHEPogdsrKmhz+vnRCovQFnaQqK5r5lgqBBzkBAf4eCAqmk5QXhg0Pkat1oXefebJeioqVWcyEeD33XBJCFAGsqmpGVlY1Ll0s5JSms9MgBAWTCPJGaCh9r7dhRZCgfUTsCIKVeHu5ssiwJ1RYWl7ehIryZlRUNnM6gqAamtAwb0RE+GDEyHAuWBUEW0BCKCzMm7dnC9VZBFW2IDu7BpcuFXF6lDrJwsJ8EB5O+6Mv13aJwBbUjogdAUZPT1kLdRy1NLfb7PXoxFFe1oTSskaUlTaxkHJ16z6R0Iln4uBArqWRE4mgBJTeCg2lqI43UhD61H5LEcay8iZkZ+VzapeijCHBXggP90FEpA+nT/W830q6S3uI2BEEK6Er2p6Opb5CqYKyskbuoiktbWQPGmrd5qvjSF+kpYVJtEbQBJTKio7x4+3J+jGqGSPxfutmGf+bLArCQr0RFe2HqCjfVxZYC4K9EbEjCFZi7ZUq1UvQ1W5xcQOKixpQXW2Cs/MgTkHRgT99RLgYEwp25bXu0jaG6oN6okC0f/cIoMrKZhb4Z07nobGxnSOXkZG+LH6iovykBkhwGCJ2BLuhxVRVf6DaBipcLiyoR0lJI4f5qVA4OtoPY8ZGITjYU9chfUF9+6OTs/L7GwkgSmvRNmZsJN9Hn42y0ka+ELh+rYTNFYNDPLmDjDa9RTdfdgyU9JbjEbEjCH28YqYiYqqxyS+oY3FDkRwq1IyL9eeDOo1AEASloLEVap14TpGcuPgA3gj67FDKq7CgAceP5aK5pR2+Pu7d4ifOj32B5EJBsAVyVBaE10Dihq5EC/Lrue5m+7Z7bMQXF+ePSZNiDT+DSFAX7W0dmhHcJGSopZ22nugPpYCLiupx5UoxamtN8PZyQ0JCAG/kIyQI/UEbnwhBcCB0tUnt349za1FY1MDVD1SMSeZ8NC6CDsp0cBYENUJDSGn6ulahFHBqahhvBHlM0cyws2e7O7/8AzxY+MTHB+gu7SXYDxE7Qp/RYy1OQ30rch/XIu9xHUytZq4zoHlQ4yfEPBW5KSpu4EJLETuCWiFBrtY0Vn8gj6m09DDeiLpaE/Ly63DsaA7bNYSFe2PI4CDExvlzu7wWkNZ1xyNiRzAk1ApeUFDHZmk0bdrXzx2DEwOxYGHSKzulyECNJoALglqh+W0eOh4NEhDoidG0jY7kKCw5iudk1+Dy5SK+MElIDETSkCCOAAnK4OzsjM7O/tl02AsRO4JhoHB4dnY1cnNqYe7o7C4oHkMpKesN0MgrpKCg3u5rFYT+YjJ1wNPTGId2+tz2dHz1CL28vFqcO5eP+vo2hIZ1z4ajlJdWoj56wM3NDe3ttjNgtQXG+EQIhkxV0VVfaSm5vFajuKSBCx2HJAVh8ZKh/R6KSWJHIjuCmjGZzDwzzYh4eLhg+PBQ3ujzX1nZgpzsai529nB3wdDkYBY/9DytH4fVnOJyc3NDW5u6jpPa+IsLgpWQkRmlpx5mVrGZX2SkD4YODca06fFPDUbsL5TGogiRIKg5suNhkMjO66I+PTO/Jk8BmpvbedDp/n0P+TiRODgQycndw3MFe0R27DtHsK/IJ0LQRf3N48e1LHAo6kIt4ePeiLJLETGFwsm0TRDUSqvJ3O/IpZ6hzi2q86GNzA3pmHH6dB6aGtvY1yclNZRnegkDx91d0liCytBquoq8b6i4+NHDKs7TU+fU9OnxUpQoGB6au2aUmp2BmBsOGxbCG0V5CgvrcfVqMWprTNzVlZoaiqAgdQsfNbszu7u7o71d0liC0C8oopKbW4sH9yu5PTwpKRhz5w3h1lRHh8fpAGmLtJgg2ONCQIwurYc+x1TATFtPGvzypSLU1bWy03NqSigCgzyVXqam8PT0REuLCWpCxI6gaqjIkJyL792rQH1DK7eHz5qdqGiencLhlP+XKc6CWpERC/0XPgkJgbyR8MnPr8OFi4Xsw0Ut7WlpYVLjYwVeXl6orq6GmhCxI6hS4NBgzXt3K9hDg8LKkybHqCas7Ofrxk6uInYENeLoied6Fj6UHqeNospU43PyRC47VA8fHsIdXzK1/eVip7CwEGpCxI4B0EpdDnU53c0oR25uDcLCfZA+Ihzh4d6qu0qluqD6+laeai4Igv6hxgRKm9NGDtWZmVXYtfMBt7Cnp4dx1EctaW01tKt7eVEaqwVqQsSOoCh0xUTtoBTFoYMFCZwJE2NUc+B4EQEBHigqalB6GYLwws/TICnXsSs0ZHXUqAjeKL11924FLl4sRHiED98XGiqjZLy8vETsCELPoM3bt8tRU92CpKHBWLwk+ZVjGtQEGbaROBMEtUGdiZ4e2vgc6QGawj5lahwmT4nl1Pu1ayWor2tFSkoot7IbNc3l6uoKs7kDakLEjo5Qe7qKvC1IJGRmViIk1AujR0f02rxrCR9yURZjQUGtYkfazh0OpdojI315o2446hjdseM+/P3cMXpMJN+vBswqSHEphXwqBLtDUZybN0pRV9+K9LQwvLk+XdOtsZRis3RZlF6GIDyHqcWs6yGgWsDV1RkjR0XwVlnZjFu3SnHq5GMkDwvhbi6tjKrQG/JbF+wCXd1kPqjilvHAQA+MHRfFtu26YVB3Ok5txdOCsWluMcPbW8SOWqD6nfnzk9jlnQxQd+9+AD9fd4x7I1pfx0OVTzwnROxoGDWmrepqTbh2vQSVFc2ct16zNlWXeWsfH3fuHpP2c0FNiP+TOqFIdmpaGG8U7aHansaGNo7+JCcHq6Ihw/zE+WQgKS1fX180NjZCbYjYEQYMRTjIbv36tRL+0L4xPhrz5g2BnqGOLCpElBOLoCZams2arIMzWrRn8eKhXF9153YZNn18h718qLZHK00ar8LPzxcNDSJ2BJ21ud6/V4mMjHKeLk6jG4ziLkodWWQnHxPrr/RSBOGpyA45fAvqh2p3yGZj/IRo5OTU4MD+Rzytfvz4aE0LVj8/P4nsCPo5oN64UYqC/Dpur3xzfRoX5RkJiuyQo6ogqAmp2dEeVPfXY1hYXd2CK1eK0dTUhjfGRSMhMUBzdYG+vr4oLS2F2hCxI1gNTQS+dKkQTc3tGDc2CtOmxWnug2graDDgjRslSi9DEJ6iw9xluAsPPREc7MUpLrqgpLIAMiscOTKcLyrJxVkraayHDx9BbYjYEV4LmWVdvtQ952TS5FhNh1htBXmZtLSYlV6GIAg6hFKRM2YmcFfrnTvl+GRTBpKSgriuhxyc1Yy/vz8aGuqhNtT9WxOULTouqOeQqrdP9wcvMNBT6WWpBqNGtARBcBwUpRs3LgpjxkTi0aMq7Nh+n2fyUZ2Pp0r9lFxcnNHRIa3nggZETnZ2DYdQaQjnosVD4eMjBY8vgg42FN3RQweFoH26uiwyF0unUJcrTVkfNiwEj3NrsXt3JsLDvLnAWQrSrUPEjvCUyLl2tRixcf5YtTpFnD5fQ1CQJ9cxidgR1AALb5Ve7Qu2iygPHhLEW35+Hfbvf8SmrZMmxarCBsPb2xtNTc1QI3I2Mzgkcqjt8eqVbpFDJoBqzwmrqZiQuieiY/yUXoogSNu5wYiPD+CtqKgeRw5nc7nBlClxitp/BAUForZWnV2qclYzKE+JnFh/rF6TKpGcPhIU7In79yuVXoYg9BoKStu58YiJ8UfMOn+Uljay6CFbDJrCroTwDQwMRE2NiB1BJeTl1XJLY2yMiJyBQAXblMYSBDXQ3NIOL4nsGBaarL7uzTRuLNm37yF3zVJ6y5HH96CgQDx6lAU1Imc5A1FW1oizZ/M5/bJqVYpqq/m1As38otZQQVADzTQqQqfDJQXroXKE9bHpbHq6c8d9xMX5Y/yEGIfMKAwKCkJNTQ3UiIgdA1BTY8LZM3lwdXPGwoVDDTPSwSHI9HNBJTQ1tmHIkCCllyGopZB5cBDP3MrKqsa2rXeRPCyEW9jtaU7o4eGB1tY2qBEROzqGpnKfO5vPA+emz4jniI5gW/z9PFBf38Z5ckFQksZGmnguaSzhadGTnBzCoyhohiGZE9KgZpq0busLNHo9uvBTKyJ2dAilVi5fKkJxcQOLnKgo6RayF6Fh3qisbBaxIyhOW3uHQ1IVgjZ9ekaNikBqaiifG27eLMX06fFsUGjL4mS1dmIRInZ0BKlqmkJ+61Ypq/epBp5d5ShCQ7xQVNyAoUODlV6KIMjnXXitI/O06fEc9af6TfJVmzkzAQE2cMcPDQ1BZaV6u1NF7GgYV5d5vf+Ojo7Cj3/ydSQmBGLDxhFwcRErVUcQEuqNW7fKlF6GIAiC1fj4uPHA0aqqZhw/nsuDjZcv+zLa2/s/7y8sLBR5eflQKyJ2NA5NmJ0zZzbMZjNWrBgupmIOhtyTTSYZCCooS1tbB9zd5HAu9I2QEG+sXZfGnmsbN27AjRs3cffuvX69VmhoKK5evQ61Ipf/GoUq6qdOnYKlS5fg4sVL2L//oAgdBVMHNJdIEJSC0hI+Upws9JMhQ4Lw0Ucfw9/fD2+/vQHh4eH97MRqhVqRSwENkpAQjxkzpuPmzVvYtGmz0ssxPDSbprbWJN1ugmI0NrSpYjaSoF06O7tw/vxF3LmTgblz56CtrQ2nTp2GydRq1cV3V1cX1IyIHQ3h4+OD+fPnwmQyYfPmrbwzCsoTGtrdkSViR1CKxqZ2+PpIZEcYOI2NTdi1aw9iY2Owbt1aZGTcxa1bt1/5PaGhYaioqICakTSWRtoGJ04cj1WrlnPK6tChIyJ0VNd+3qL0MgSDGwpKGkuwJYWFRZza8vT0wFtvbUBw8MsNKyMjI1Baqu5GDRE7KiciIhzvvPMWV8l/9NEmlJWVK70k4RlCQrxQVdms9DIEGD2yI2kswbZ0dVlw8eJlHDx4iFNbM2dOh7Pz815OUVGRKC0thZoRsaNSXFycMXfubEyZMhk7duzm+hwVm1MaGnd3F7S1y4wsQTmaGqVAWbAfdXX12LJlG6qqqvniOy4u9qnH/f39+TlqRmp2VAjtSLNmzcTly1dw/PhJpZcjWAH5GpFzNZl2CYKjMXfIvifYn3v37iMnJ4ejPCNHjsDRo8dhNrerekxEDyJ2VIS7uxvvRE5OTtiyZatqB6oJzxMW5o2Kimab2q8LgiCojdbWNrY6oYvyjRvfREbGPVU7J/cgYkclDBkyGNOmTcGZM+fw+HGe0ssR+khEhA/KyppE7AgOh66qZUyE4GgKCgrx8cebsXHjev7azc0N7e3tUCsidhTG1dUV8+bN4X/TjkNOyIL2iAj3QXZWgdLLEIxqKCiGooICmM1m1NXVIScnF2+9tV7VF+sidhSE5lmR0Dl37gLvLIJ28fVzR0OjpB0Fx9NQ3wY/f+nEEpTBz88P9+8/QHZ2NubNm4uUlOE4duz4gOZs2QMROwpANTnUwhcQEIDNm7ep2mJbsA5KIzh9OjaCfJEEwVHUN7TC389D6WUIBsTDw733/EXi5sCBQ0hMTGBfHipeLilRTzu6iB0HExISgsWLF+DGjVs4efK00ssRbO23U9XCxcqC4MjITmSUr9LLEAxITEwMioqKnrqP0lhlZWVYsmQxysvLeQSFGrq1xGfHgYwbN5bTVrt27eUWPkFfhEf4oLysSellCAajvqEN/n6SxhIcT2xsLDstPwvN09q+fSeam1u4gNnPT3kxLmLHQaG+NWtW8e3mzVvQ2Nhol/dxdZnXuwlKdWTZ528rCK8cAipiRxjA+cK1n+cMcvh/las/meEeOXIMK1cu51oeJRGx44Ai5A0b3sSVK1c/DecpvSLBXgQFeaKm1qT0MgSDQXViNHVaEByJ06e1ia+bdl5dXY2PP/4E0dHRWLZsCVxdlamekU+InSDbi6lTJ2PSpIk8obyoqFjpJQl2hr1OLN2+J4IgCHomIiLC6lmNnZ1d3KH14EHma4eK2gsRO3aApsSuX/8mTyanvKU4IRuH4ODuImVBcATt7TQmQg7jguNJTEzE48eP+/Q9ZLGya9ceLFq0EKmpKXAk8imxMTTqnoTOmTNnce3aDaWXIziYqGhflJRI3Y7gGBoaxGNHUIbY2Oc7sayhoaERmzZtRkxMNBYtWuCwFKyIHRsyevQozJgxnafDlpaWKb0cQQGiovxQUtyg9DIEg1Bf3wo/8dgRHIyzszOn7Ts6Ovv1/VTnQ4XLeXn5eOutjfD3t/+YHfHZsdEffvHihWhqamKhIzUbxsXf3x319ZK2FBznsUP7nCA4kpiY6H5FdZ4lM/MhyssruFvr9OmzyM8vUFbsrFmzps8v/Lvf/Q5hYWEwglX2ihXLcPnyFWRlZSu9HEFh6GrHxcUJZjPVUjgrvRzBAO7JYigoOJrExARkZ+fY5LVqa2vxySdbsHz5MtYMV69eg2JprF27dvFEU39/f6u2/fv3c5TDCOp21aoVPO5ehI7w7AR0QbA3EtkRlCAqKhIlJSU2ez0aNUHNPB4eHli+fClnSxRLY/3v//2/rY7UbNu2DXpn5MgRGD58GLeVU9eVIDxVpFzcgNhYf6WXIugcU6sZnp6uSi9DMBCurq7cSk7+Trbm7NlzSE4eyq7Lu3fvtWnQxKrIzsmTJxEUZH1f/MGDB9lASI+QlQqNfCDnyK1bt4vQEV5cpCwdWYIjkPJAQYEUVm5u31rO+8KjR1k4fPgITx2g7maHip2ZM2fCxcX6WuZp06bB3V1/oVVK5a1btxZVVdVcSS6FyMKL8PBwQVtb/7oUBMFapC5MUIKhQ5OQnW3fsg06x27dug2zZ8/kSI8i3VgNDQ0vLcwkgUOCQI/4+vpyfc7p02dQUFCo9HIElePt7Yqmpnb4+Ojz8yAoT11dKwICpe1ccCyBgQGora2z+/vQMFEqE1m2bCkCAgJ45JJDfXboTQMDA5/b6H5PT0/Ex8fj/ffff+28DC0RGhqKNWtWciGyCB3BGmJi/FFUJH47gv2orTUhIEDEjuDYc2FFRaXD3o9qg6h2h6YSkAEhj+RxVGTnj3/8I77zne/gs5/9LCZMmMD3XblyBX/605/wT//0T6isrMTPfvYzjvJ8+9vfhtZJSIjHtGlTsG3bDh5XrxWenGJr7jim6FqMSGycH27dLMPw4SFKL0XQcWQnNMRb6WUIGsK1n9PNn0xhKdF5TB481BS0bt0a7N69h7u37C52SNT8/Oc/x/r163vvW758OUaMGIH/+I//wPHjxxEXF4cf/vCHmhc76elpPJaeQmlmc4fSyxG0NiOrWjviWNAedbWtGDo0WOllCAYiPj6OPeWU4M6dDNTX12P9+nXYsWMXWlpM9k1jXbhwAWPGjHnufrrv4sWLvQXKBQX2c0J0BFOmTEZcXCy2b98hQkfoMxRudXVx4kGNgmAvQ0E/P/01ggjqxN3dHR0dHejsVO6YRg7LR48e50ahvo6Y6LPYiY2NxR/+8Ifn7qf76DGiurqa63hsza9//WskJCSw8dDEiRM5fdbDw4cPMXXqVMTExOAHP/jBgN5n3ry5cHFxxoEDh+ziJSAYg6hoakGXuh3BPtCxyVFDFAUhKWmIzVyTBwKNl9i7dx83DAUFWR/Z7PMnhepxfvnLX2LUqFH4/Oc/z9vo0aPxr//6r5zeIq5evYoNGzbAlmzevBlf+9rXuPj5xo0b/P4LFy5ERUUFP/6Vr3wFn/nMZ7B7927eKALVn6vxpUsXo6GhHmfOnLPp+gXjERvrh8JCETuC7RHbC8HRDB8+DA8fPoIaoG6wbdt2Yv78ufYTOytWrEBmZiYWL16Mmpoa3ujfdN+yZcv4OV/60pfwi1/8AraEXu8LX/gC3nvvPaSmpvLsLS8vL3zwwQe98zXGjRuHkSNHIioqCnV1fWuNGzTIiZViUVExrlyxz2wOwVhERPiitFTMBQXb021rICkswTG4uDizrUxLi3rqEJubm7Fz5277Tj1PTEzET37yEziK9vZ2XL9+Hd/61rd673NycsK8efN664S+//3v89cmk4lFF0V9+gJNXX3wIJOnsAqCLaCBoJYuC6cbnJz63zIpCC8qTg6UtnPBQQwePBg5OblQG+3t1k8w6FfC9+zZs5wymjJlCoqLi/m+P//5zzh3zj6pn6qqKi6KCg8Pf+p++rqsrIz/vWTJEm57p+FkO3fu7PMgsVu3bovQEWxOWLgPystlKKhgW2rrTGIoKDiMlJThmj8/9lnsbN++naMmZCBItTM9s6GoJexHP/oRlK4WJ9OjvtDTHp+Xl2enVQlGr9spKqxXehmCDiM7AQGeSi9DMABOTk7w8fF56fQE3Yod6nSiepnf//73PP20B+qEIvFjD0JCQjhSU15e/tT99HVERP8Hhf3N3/wNd3gZwUjqyU1wrJNyoTgpC3ZwTw6UyI5gBQM99sfHxyE/Px9ap89ih1q8Z8yY8dz9/v7+fS4KthYqjKLiYzIs7IHGUdDXkydP7tdrkgP0b3/7W+7uEgT7DgXtkO4Zwaa0tnbA0/MvF5uCYC9SU1O4ntVwYociKS+aeEr1OlTEZC+o7ZyiSeTg/ODBA+74omps6s7qK4cPH+bXoXZ5el1BsCdhYd6oqGhWehmCThDhLDgyhRUYGIDq6hponT53Y1H793//7/+dW77Jl4YKgqkj6utf/zq++93v2meVAPv2UAHy9773PS5KJm+fQ4cOPVe0bA1Uc0St8sOGDdN8HlJQP/HxAcjPr0N4uI/SSxF0QEuLGV7ebkovQzAAgwcnIifnMfRAn8XON7/5TU4hzZ07l3vuKaVFhcEkdr761a/CnpBxIG22gISOIDiC2Fh/3L5VhgkTYpReiqADaqpNCA6S4mTBMfMhT5w4CUOKHYrm0NTzf/iHf+B0VlNTE5v8UbW2IAjP4+7uAnNHl/jtCDahuqYFQcEidgT74urqwl3XDQ36MEbtl6lgT9EwiRxBEF5PRLgPysqaEBXlq/RSBI1TU2NCenrf0/eC0BeSk5Px6JE6xkM4TOysWbPG6hfcsWPHQNYjCLokISEA+Xm1InaEAVNbI23ngmO6sPbtOwC9YFU3FrWV92x+fn7c8n3t2l/mR9EoB7qPHhcE4XmiY/xQJH47gg3o6OiCq2vfHOIFoS94enpwyQqNXzJUZOfDDz/s/fc3vvENrF+/no0Fe0Yy0CgHMugjISQIwvPQyanLYkFnZxecnfs1pUUQuO6LTkKCYE/S09Nx9+496Ik+1+xQyzl56jw5e4r+TX41NCvrX/7lX2y9RsHGPOmkae44puhajERUlB9KShq5O0sY+EmfIhwkHvm2owsdnZber2kAa68bjYX++/SrT28GOQ2Cs/Mg9hGhWxKgPbdURO7q6sSDXNUmLBoa2uDnL9POhZdjC5f8YcOSsWnTZhha7HR0dPR61DwJ3Uct6VqF2ufb2tqVXoagYxITAvD4ca2InU+N8czmLrQ0t6PFZGbvmJ7NRJupA+3tnby9iEFOgItztyBhoeLiBBcSK3zbLVi6n8j/URtpz5f8PxJDnZ+Ko07+NwmnT+8j4fTp9sL3HtTdYefpRd0qrrx58S3d1/21t7erXSJ4NdUt0nYu2JXw8PDe4duGFjvkWPy5z30OOTk5mDBhAt93+fJl/OQnP+mXm7FaGDp0qO7CdoK6iIr2w7nzBdA7JBwaG9rQ2NiOxsY2NDTSv7u/7jD/RUC4uTmzOPDydGGTPBIMERHu8PJy4zEbtFGERW3RFYoq0QgQ06cijYQZCbSGhla0fPrv5pZ2fh4xCIPg5e0KXx83+Pi6w8fHDf7+HggI8OCfsS9U15gQFOxlp59MEIAxY0bh5s3b0Bt9Fjs/+9nPeGTEz3/+c5SWlvJ9kZGR7Lvz93//99AqyckidgT7QtEGdzdnPklqfa4RneTr6kw8fbu2rhV1tSY0Nbf3/px+vu7w9XOHr687IiJ8MHRoMP+bBI7WoZ+vJ6ITFPT655Pood9XU1O34GtqbEdpSSPq6lrR1t7Bz6FoFAkg/4BuERQc7MlTzZ/1ZaK286QkK95UEPqBs7MTD95+dui2IcUO5bj/8R//kbeeUQt6KEymFBxVoJtMrUovRdAxCQmByMurQ0pKKLSQamqob0NVdQuqqro3itAQFIWh9ueAQE/Ex/lj1KgITt2oLQqjBkiwUDSHtoiIFz+HUmb19a0sgGijdCfdUrqNUnOUuqKITllpI6fQBMEeJOvMW+dJBvSp0YPI6eHRoywMHz5Ml+E7QT0MHhKIc2fzVSd2aIp2RUUTysubUV7WxFEa0i3+fh4IDvFCWKg3r9nX100EjR2g2qPgYC/ensVs7uSIDo2JoAjRkcPZaG3r4GhQSIgXz1wLC/dGYODzkSBB6AsjRqRj7959MKzYGTt2LPvoBAYGWvWi06ZNw+bNmxEdHQ2tkJWVjbVrV4vYEewKpSqoo4aiJkqJBjphFhc3cJSgorKZC4UpvUYnzfAIHxY1FIUQ1AHZFtDfhtJbmZmVWLU6pTcaVFnZjIqKZty4UcKCyNJFF6HubF4ZFe2LkBBvEUCCVQQE+MNsNus2u2GV2Ll16xZu376NIGsS1J8+v62tO9ytFczmdu408/b2QnNzC4yCtKE7HoqUVFe38InI3tAJkSI1xSUNKClu5AiOl5crIqN8kTg4EBMnxeqijsYI0D5D+86T0aDISF/enko9NrSxxcHdjApUVjVzuz19X1SkL4sgqguS6Jy+sEW7+bhxY3Ht2nXoFavTWDTlnD5I1qDVDxIVKJOZ0uXLV5ReiqBjhgwJQk5OrV3EDhU/FxbWoyC/nmtsyDuGojXR0X4YMSJc84XRRqaqsgWhr9ln6NjLhc7+Hr2pUjpuk1AisXvpUhFq60zwcHdBTKw/2yCEhUn0x+g4OzsjKioKx4/rY8J5v8XO48eP+/zCMTEx0BpZWVl4662NInYEuxIX589ph4kTB/4ZaW5uR0EBiZs6TmOQmImlguHREVzPodULD+F5SLyOGNn3AaC0D5Cwpm3kqIheUUzjS+7fr8CpU80sisn0Mj7en29F/BhvDtb9+w+gZ6wSO/Hx8TAC5A9SXV2NsLAwVFRUKL0cQafwXCNLd+FpX2cckcleQUEdcnNrUV3VAm9vNxZP4ydEc4GqiBv9Ul3TgiAbGQqSKCY7ANoI2hcp9UX71blzBXB1cUJcfAAGDw7k95T9St+MHJmOrVu3Q89ID+Mz3Lp1G6NHj8SRI1K/ItgPir4UFtRj8JBX18FRCqK8vIlPQkWFDXzFHRfvj7FjI7lzR05CxqGr02K3uWokuuPjA3gjyDQxP78O16+VsMjy8XFHYiKJnyCu+RL0Q3h4GGpqatHeboaeEbHzDKWlZZg3by6fVHocUAXB1pAx3I0bpS8UOxS9yc2tQdajajQ1tXPNzeDEQEyYEMNFqYLxoMiLi6vj/vbk5ZOcHMIbQUXPeY9re9ve42L9MTQ5WFKlOuCNN8bhypWr0Dsidl5ATk4uBg8ejOzsHKWXIugUispQGqqnBZ3M+kjcUASHhlbSFfTMWYncRiwI3Iml4JgI2g+p3oc2SvdTVPLWzVJUVnZ3iA1NCuK0l4hxbeHh4c5+eZWVVdA7gyzWtljpFHKB9venwYz0Ie2+QvHx8cHChfOxfftOGBVpQ7c/e/Zkwt3dmc3iqPYmOTmY28HFIVd4lrsZ5SyK09LDoCa6O71MyM6qRl5+HTtrD08JYbEuwkf97eZTpkzmOtWHD7XqmkzyhdzH619rcixH1RfQ1NTEBxYfH280NTUrvRxBR9TWmvDgfiXXQzjRlG6nQdiwcYR0vwivpLKqRXWu23/p9PLibdLkWE53PcysxPZt93jI6fCUUC5y7mshvmB/nJyckJQ0BBcvXoQRsJnYeffdd1FYWIgTJ05AD9y8eQtjxozG2bPnlV6KoHHq61rx4EElzzvy8/dAakoonxjoqnjb1nsidITXUsNpLNt0Ytk73TV+QgxvlJrNzKzCjh332dcnNS2MhY+9iqyFvpGWlooHDx7AKLkdm4kdGg1BSlEv5ObmYurUKTh37oLVZoqC0AN1s2Q+qELmw0oO7dOBntrDnz7QD+KUFY1vkA4X4XVO2FqLjtCU+/Hjo3kj4XPvXgWuXS3m2qP0EeGIjPSR4mYFGTVqJDZv3gKjYDOx86Mf/Qh6gvQNFSoPGSKFyoJ1UPdefl4dMjLKuWMlZXgoVq9OfeU4BnZTzq7pl1mcYByho/VoCAmfSZNieaOBs/QZOXXyMeITApCeHsaOz4LjiI+PQ0lJCczmDhgFqdl5zYyvJUsWi9gRXkltjQm3bpWitLSRfUpmzkqw+uA9JCkIhw9lidgRXkoNmQlqIIVlLWFhPpg714cvDvLyanH6VB7a2jt4nAmZHGpd2GmBCRPG4+DBwzASfRY7a9euxYQJE/CNb3zjqft/+tOf4urVq9i6dSv0Ag0E7ezs5Cpv6toShB7oQP3oUTUy7pTBw9MVY8ZEYNbsxD6H5Sl9Rb461M4rB3nhRdBUc5pfpTeoVo26tmij8RV371Zg8ycZPNh09JhIdgQXbE9wcDDa29u5EcdI9Ln1PDQ0lIuQR4wY8dT9GRkZmDdvHsrLy6H11vMniYuLxdChSboekPYipPX8xTTUt+LmrTIUFdazqZothmtevFCAiIjuKeSC8CwnTuQiPT1cl4LnWeh0VFzcgJs3y9BqMnNtD1kyyIWA7VrPly9fyvMfKyoqoX3s2HpOatDNze25+11dXXUZ/SgoKMSMGdPh6upiqPym8PQBOC+vjod3Ojs5YcyYSMyYEW+z4koSTTeul4rYEV4ImU9qoRPLFtBnKibGn7fW1g6u7flkUwYSEgP5cyeF/APD398P7u7uOhE6faPPcpkiOps3b37u/k8++QSpqanQI3fuZGDkyJFKL0NQoDD09u0ybPr4DjvGLlw4FKtWp3BRpS27SKg7hSZaS9ef8CI6u+w3E0vNkE8PdXK9/c5IhId7Y//+hziw/xEqK8X7bCAmghcuGMNXZ8CRne9+97tYs2YNcnJyMGfOHL7v+PHj2LRpk67qdZ7k3r17eOedt3Hjxg3DeBI8GyI1UlqLWsFvXC9h4z9qGX9zfbpd235JOFEbbmlpE6KifO32PoI290VPD2P3kdDnIykpmDcSOlevFvPMOIr00Iw5I7WvD8Q12dvbm0s2SkpKYUT6/Clavnw5du3axa3m27Ztg6enJ0c9jh07hpkzZ0KPUPFoXl4ehgwZIp1ZOu96uXSpCM3N7Rg7NgpTp8U57EA6bHgoO8+K2BGepJKKk8N9lF6GaggN9caSJcksAm/eLMWVy0XcyZiWFmbI6FdfmDx5kmGjOkS/LhmWLl3Km5G4du0GF3aJ2NEf5Ptx8UIhCxtyNlaiEJQiO6dO5vYOBhWEnn3TCIXJfYVqd6ZOjcPEiTGcaqa6nuRhIRg9OkJz5ouOwNPTA2FhITh27DiMSr/ETl1dHUd1yGX461//OoKCgjjFEx4ezk7KeqSlpYW30NAQQ0yINQIlJQ24dLEI7h4umD4jAUFBnsrOGAr15todunoVBKK8opnnSwkvhoaNjhsXxSmt+/crsHXLXfa6GvdGNNf8CN1MmDABly9fhZHp895w584dbjGn3B+ldj7/+c+z2NmxYwcKCgrwn//5n9ArFy9e4gKvPXv2Kb0UYQBQ2/jFi4Vs/Ddv3mCeV6UGhg8PQeaDShE7Qi/NTe3w8Xm++1V43rOH2vMpnZWdXYOdO+4jNs6fC5xpJIuRoe6r2NgYnD59Bkamz0nOr33ta/jsZz+LrKwseHj85SSxZMkSnDmj719mVVU1t9h3+/IIWoMcjmnwJg3lXLwkGQsWJqlG6BDUbltUpD/7BqH/xpUU8ZO0pvXQ74pcmDe+NYLTfzR9nXysyLjTqEyePBGXLl2G0emz5CWX5P/4j/947n5KX5WVlUHvUIHXlCmTDGe1/WQXgNY6s6iD49zZfE5XscDxc4dar079Azx4/ESggik1QR3U1ZkQGKgeMa410ZOcHMLC52FmFae3koYGY+zYSE3V9Ayk+4qgBqLo6CicOqXvQIRdIjsUEnuReeCjR4/YXVnvlJaWsVOjt7eX0ksRXkNdXSv27snEpYuFmDUrkbs41Cp0ehg+LASZD6UmTAAqypsRKsXJAxY9VPNEXj2UDtyy+S7u3CkzjKfV1KmTcf68cTuwBiR2VqxYge9///swm829OxPV6tCsLJqbZQQoJDhp0kSllyG8BHJePX48B8eP5XB31fIVwzUTKSHDQpqcLghUnExDM4WBQ+cpqueh9FZbWycbhebm1kDP+Ph4cwAiLy9f6aVoU+z8/Oc/55ERYWFhMJlM7K2TlJQEX19f/PCHP4QRyM+nWUbh8PBQd5TAaJAfEhmO7dh+j4cLrl2XprliX/IK8fZ25RlcgrGpqmxGaKhEkG39+aKi5TVr0/iiYuvWuygv1+dAzGnTpuLs2XNKL0O7NTtUnHv06FGcO3eOO7NI+IwdO5Y7tIzElSvXMH78eNmZVACFpLOyqnH1SjF3ZGx8ayTXv2iVlNRQLqKeOClW6aUICmLu6NJUfYmWoLb02XMG80XF6TP5fLyYNSsB3t766HyjUgvaioqKlV6Kauh3T960adN4MypZWdmYNGkCrlxxQ1tbu9LLMSzV1S08FTo8zIfHOri5af/kkJgYyBEqETvGxWSSMRGOgLoxly8fxpPW9+zO5PET5NGj5YslYvp0ieoMKI3V1dWFDz74AMuWLUN6ejoPBaUaHvLWMUrB15NcunRFancUglpJT558jNOn8zB/3hDMmJmgC6HTE2r39/NATY1J6aUIClFe1oTwCKnXcRTR0X5cz+Pq5sz1PDQXT6uEhATDzc2Nm2mEfogdEjMkbMhEsLi4mIVOWloa8vPz2Xdn9erVMGJ0JyYmhq24jQS1Q/Zsjob2w8zMKu6qiInxw5o1qQgI1EbxcV9ITQtlR1jBmJSVNSFCxI7Di5hHj47kWr+HD6uwe9cDNDa2ae74Onv2LJw8edrm69I6VsdJ//jHP7JpIE04nz179lOPnThxAqtWreIIz1/91V/BSJDvzuTJk3HixEmll6J76mpNOHacUlbefBVGVvF6JS4uABfOFwLGzRTD6GJn9JhIpZdh2HqeBQuSeC7Zvn0PMWxYCI+j0IK5Y3x8HOrr63mkk/A0Vp8tNm3ahG9/+9vPCR1izpw5+OY3v4mPPvoIRuPx4zyEh4eJ746dnWSvXCnC4SPZmDNnMM+x0rPQIahmICTECxUVzUovRVCA1rYOme2kMNT2v2HDCHR2WrB5812eW6dmSItNnz5NanVegtVnDOq8WrRo0UsfX7x4MW7fvg0jcv78BUydOgVG5MmQqz3SWuR+vPmTDO5KWb8+XdFhnY5GUlnGpK2tQzf1Z3q46KBW9SWLh+L06cc4eyaPLS5sjS2OoVRWkpOTC5NJbCsGJHZqamp4qvnLoMdqa2thRAoKChEYGMheQ4JtoAPKmTN5OHsmH0uXDdNMGNnWRZMlxQ2GLP43unNyeLjU66ita4vqA0NCvPHJpgyes6cmnJ2dMW7cGFy5YuzJ5jYRO52dnXBxcXnlL7ujowNG5fTps5g1a4bSy9AFlCunaE5oiDdWr0lR/YgHe0HiLiLCl+s3BONQVi7FyWr9PJIH1uo1qbh8qYgvxuwR5ekPEyaMx/XrN/k8LbwYq5PCdHVJXVc0G+tFtLU5vmpdTdAQVCcnJ3aWrqiQ1EN/oH3s8uUi9rxYsTKFZ9kYHUpl3btbgchIiRoahbLSRqSm6n/OoFbx8nLFylXDce9eBTZ/chcLFg7hiI9SUL3okCGJ+OijTYqtQVdi5913333tc4zWifUsp06dxqJFC7F581all6I56utacehwFpKHhnC42Ggpq5dB6YwTx3O5SFvrRmeCdbS0mHXj5KtX6PhEbu3UNXn4cBbfTpgQrchxa86c2Thx4hQk220jsfPhhx9a+1TDUl/fgKqqagwenIjc3MdKL0cz0Zy7GRW4e7ccCxcNNVQBsjXQwTM+IZDn+CQODlR6OYKd6ejogrPOOw31BKXY161Lw40bpdi69R6WLEl2aEQ6KqrbnqCkpNRh76lV5FNlY86dO8+dWRKZsM4Fed/eh6irM2HDxhEidF5CenoYi0HBGPVqYRobXmt06Fg/blwUZs9OZCPC3BzHTFOnUwwZCIrHm3WIkYONodqlhw8fYuTIEbh9+w6MxpOtk+aOY688qB89koNp0+MRHx/goNVpE39/D5hMHSwOpSVZ35SVNiEiUoqTtUhoqDfWb0jHsaM5yMuvw8yZCTz65UXYwqZj5MiRyM7OQXOzuv1/1IJEduzAtWs3MHr0SLi5uSq9FFWmrW7dLMXpU3lYtTpFhI6VDB8egszMSqWXIThiTIS0nWsW8gNbvCSZu+lopE1dnX08b9zd3TBq1AhcvXrNLq+vR0Ts2AEamHr+/EVMnz5d6aWoMm3V1NyOdW+mSRFmHxhGYudBldLLEOxMQ0MbfA1qtaAnUlPDsGjxUBzY/wg5dkhrzZgxHefOXeBzjWAdInbsBIUXg4ICERQUpPRSVAFN8N665S5GjozAtGnxUtPUR9zdXXh8QEO9uKPqFbO5Ey6uTvLZ0AmBgZ54c30aHtyvxPlz+TYzByV7Ez8/P2mC6SMiduzIsWMnMH/+XBgdKtg7dDALy5cPQ3yCpK0GVqgsHk56pZzMBCWFpbu01tJlyXD3cMGunQ94FMhAIB1M55QjR47abI1GQcSOHaHxGRUVlUhOHgojQlcyly4W4t79Sr7CIct1of8kJAYiL69WxkfolJKSRkRFiXmk3qBI3RtvRPOMrW1b76G6uv8FxWPHjkFWVhYaG8VVva9IN5adOXfuHN5+eyNyc3PR0WEcK29XVxdsWP/vKC0txaVLV17ZmSVYB5kKhkf4cBGrOCrrj9KSRowY8fL5g4K2SUxcCx8fH6xcuRxDhgzmoZ19wdvbGykpKeKU3E8ksmNnzOYOXL16HVOmTIZR8PHxxoYN67n1noSOYDtGjIjAndtlSi9DsANkL+DpKR2ceqapqQmbN2/BqFEjMW7c2D5974IF87g0QiK7/UPEjgO4f/8BoqKieDK63gkJCcHatatx6NBhKaCzA2Fh3lzsTU67gn4wmcxcgC7oH4rw79y5i4uMScBYU5BOkaDm5maewSj0DxE7DuLw4SNYuHA+9ExCQjwWLVqAbdt28tgMwT4MGxaCh5nShq67ep1oSU0aBQrOnDx5CuXlFXxx6Or68oiei4sLpk2bwrMXhf4jYsdB1NbWobCwiJ2V9Qj9XOPHv8FDUOkKRLAfaelhuHdfurL0hBQnGxNK9V+7dh0bNrzJ08tfxOzZM3HhwiW0t5sdvj49IWLHgVy8eImdlb289DUDatq0qTyQbtu2HTCb5QPpCM8dqu2orTUpvRTBRpSXNfGEe8F45OXl4/Dho1i3bg38/f2feiw6OgpeXl7IyspWbH16QcSOAyG3y+PHT2L+/IHPRVEL5PlAKedDh45I4ZwDGT06ArdvSf5eD9DnhmqwyJNFMCaVlZXYtWsvVq1ajtDQUL7PxcUZc+fOEU8dGyFix8EUF5fwsNDExARoGSqqW7ZsCerq6nD27PnXPp8G3z25CQMjJsYPJSUN6OoSgal1Gurb4O8vIyL0SF+OefX19RwdX7RoPmJiojFz5gxcuXIVJpO4ptsCETsKcOLEKcyYMe2VRWlqhib5rlq1AgUFhdxWLygjNgcPDrLL3B1BiXodP6WXIagAmmC+Zcs2zJs3lzt4MzMfKr0k3SBiRwHa29t5iNu8eXOgNSi0unbtGm6nv3MnQ+nlGJoRI8ORkVGu9DKEAUIROunEEnro6OiAxdKFzs5ObjkXbIMYOygEuWempAzndm0qUNOS0KHugb66fz7Lk2FdcVfuHzQ1nhw6mpra4eMjE+S1SlVVC4KDX9yJI2gLW6Top0+fjmvXbuDhw4d8vKUoLg2WFgaGRHYU5MiRY5yXdXNTfzrL2dkZa9asxvXrNwYsdAQbR3fuSKGyVuns7OKTGY0CEQTqvgoKCsS9e/fZfJBqeGgellHnK9oSETsKp7POnDnL+Vm1Cx0yvrpx46ZcYagMqtvJzZXhoFqedE7zzgSBajiptOHAgUO991Eqa/v2nTxeQgTPwBCxozCPH+fx7eDBiVBrMfKaNatE6KgUiggkJATgcW6t0ksR+kFRUQN31gkCjY44c+YcWluf7r4iwbNjx06MGTNa8128SiJiRwUcO3Yc06dPhbu7utpPyT9nxYrlXIgsQke9jBodiVviuaNJiosaEB0tYsfoDB2axIXJPRe/L0p3kuCZPHkit6ULfUfEjgogG3BqR1+8eCHUxKJFC3mY58OHj5ReivAKqDjZ2XkQ6uvFj0NrtLV1ygBQg0MOySRijh8/8crnmc0dnNKaNWsmwsPDHbY+vSBiRyXQ3Cyan6WW2Vk0j4XWQ7NbBPUzZmwkbt4oVXoZQh9oqG+Fn5+6ormC41myZBEOHz7GBcmvo62tnQXPwoXzEBwc5JD16QUROyri7NmzGDEiDYGBAYquY+LECdwhcunSZYe8nzgrD5zYWH/2a6Fwt6ANpF5HHwzk+EWdVqWlpSgvt94vy2QyYefO3exg/7LhocLziNhREWT9v3//QSxZslixVlTy/gkNDeG0mqAdSJwOGx6Kh5lVSi9F6IvYiRWxY1RoBhZ1WF24cLHP39vY2IQDBw5j9eqVmnXidzQidlRGXV09FwTPmDFDEY8HSqM92fooaIf09DBxVNYQNTUmBAZ6Kr0MQQFcXV24RnPfvv3or2sEDQ8lJ34a3UMXO8KrEbGjQjIy7iIgwB9xcbEOe096v7lzZ2PXrj08nV0pZGBo/3F3d4GfvwcqKpqVXorwGtrbO+Hi6iQnKY0y0GPUwoULcO7ceTQ1DeyzSu77ND9r0aIFA3odIyBiR6UcPHiIi4SpUt/eeHi4Y8WKZdi9ey9PZBe0y7hxkbhxvUTpZQivobSUhn/KPCwjkp6exnU31Olqq4vjpqYmTJw43iavp1dE7KgUqro/dOgoixB7XvxRbdCqVStx7NgJ1Nc32O+NBIcQFubDLehtbR1KL0V4BVKcbEwCAwPZDfnkydM2fd2zZ8/zlHQxHXw5InZUDFXoP3r0CNOnT7Pbe8yZM4fnsJSUSNuyXhg1KgK3b4vJoJopLWlEZKREdowEudEvXboY+/YdsEupANX/TJs2FQEBynbzqhUROyrnxo1bfDVgD8WelpbKH0AKg6oVqd/pO8nDQpD1qJq7+wT1QXPMzB2dcHV1VnopggNrCRcsmI9r166jvr4e9oBMB/fs2YcVK5ZqYri0oxGxowGoO4qiOz4+thsYGBYWilGjRuDo0eM2e01BHVBqckhSELKyqpVeivACKitbEBbqrfQyBAcyYkQ6zGYzFxPbExJSp0+fxfLly+z6PlpExI4GoA8JFSyTYreF/w4VJFPbIxUkK9l5JdiP0TwvS1KTaqSgoA6xcf5KL0NwEGFhYVyUfOLESYe8X35+AYqLS6Rg+RlE7GiEysoq3LlzF/PnDzyds3TpEi5Ibm5uscnaBPVB85b8/T1QVtak9FKEZygqbGDHa0H/0HDnxYsXYM+evQ5NK5P7fVxcHKKiIh32nmpHxI6GuHv3HkdiBjI/a9y4sSgrK2PlL+ib8eOjcfVqsdLLEJ6p1zG1muHpKTUVRmD58qU4fvykIheWVLA8f/5cjuQLInY0B0VkUlNTEBkZ0a9w6tChSf2yJxe0R3CwF7egNzW1K70U4VNqa1sRJK7JhmDq1CmcUioqUuaCw2Rq5fPFsmVLFXl/tSFiR4NXhrt378GCBfP6NASO7MnJZXPv3v7bkyuNdGb1nXFjo3DjhkTx1FSvExcvrcF6P94MHpzIMwavXr0GJaEIfnFxMd54YxyMjogdDUKK/eDBI1i5cgWcnJysbnske/LmZhklYCQSEgNQVFiPjg4pRFcDhQX1Uq+jc4KCAjF16mTs338AaoDqd4YOTUJQUBCMjIgdjVJRUYGbN29h4cL5r33ukCGD0dnZaTN7ckE70OyllNQw3L9XofRSBABNze3w8XFTehmCnXB3d+O0EXW6ku+NGqBI/oEDB7kD1xbdvFpFxI6GefAgE62tra8MUVJx2rRpU3D8+AnoCUlpWc+IEeE8DZ1SoIJy0BgPPz8PpZch2Om4QhcWFG2nFvOGhkaoifr6Bm5woToioyJiR+PQjJXY2BgkJQ15afrqxIlTqrnKEByPi4sTEhICkJNTo/RSYPQUVpz46+iWuXNn4+HDR4oVJL+O27fvIDw8HBERfW9u0QMidnTA3r37MHHiBO62epLk5KEc+SksLFJsbYI6GDsuCtevSaGykhSI2NG1QzLVT5KgULsb/4IFc62u9dQTxvuJdUhHRyd27tzNOVkfH+9eM6vJkydyVEcQyNclKNgTxcUy2V4p6hta2ehR0Bdk3Ed2IFoYvdPS0oLbtzMwadJEGA0ROzqBduL9+w9i1aqV3GY+e/ZMnpHS0aH/9JUthvQZgYkTYnDlskT5lKC5uR3eXlKYrEYGcuzw8/Nj4z6yA9FKTdzt23cQHx9nuOnoInZ0RFVVFc6fv4D169/kyE5eXr7SSxJUhJ+/B9fvVFfLmBBHIyks/UHH2FWrlvOk8dbWNmiJw4ePWtXJqydE7OgMEjj+/n6GiOgIfWfipBhcluiOw8nLq+MicUEfUM3L6tUruUygtrYOWqOmpgYlJSU8oNQoiNjRGTT76sqVq2hoaGBjKyMiKa2XExbmg5YWs4yQcDB1dSYEyJgI3aS9ly5dzOkgtXZeWcOFCxcxduwYuLkZY06biB0d4eXlieHDh+H69Rs4e/Y855NHjRqp9LIElSEDQh1LQ0MbfH1lGKNemD59KpcMkM+Zluns7MLFi5cwdepUGAEROzpixozpOH36TO/sq0OHjrB7MlmFC0IPVDtSVtbIQ0IF+5OXVyspLB21mPv4+ODixcvQA1lZ2QgLC4W/v/7ryUTs6ITg4CB4e3s/5anTPTR0Lzssx8REK7o+QT2Q0+vYsVG4fl18dxxBfl4d4mX4p+ahDqaUlOE4fPgI9MSJEyfZEFHviNjRCbNnz3qhpw7NxNqxYyc/HhISosjaBPWRnByMvMe1aG/vVHopuoYuOKg+StJY2oachyl9tWvXbnR1aaPF3FoqK6tgMpkQFxcLPSNiRydXHI2NjaitrX3h421t7dixYxeWLFlkOG8F4eXRndFjInHrZqnSS9E1tbUmBAVJYbKWCQwM4Dbt7dt3ob3dDD1y6tQZTJ8+DXrGRekFCAOHrji2b9/5yuc0Nzdj1649WL16Bbstq21QneB4hg8PxaaP7/AoCfLfUTMUgaJC36amNphaOtBiMnNXmaml+7bd3Ak8c8FtgQWD8PyU5+fup39aAGfnQew07enpAg+69XCBp5crR2X8/Nzh7e3KIrEv5D2uQ0JiYL9/bkFZqDRgxYpl2LlzD0c/9IrJZOISCBox9OhRFvSIiB2NQztnQUEhTKbW1z6X2tH37NnP/hDbtu1kAaR3nmwtNXccU3QtasPJaRBGjozA7dtlGDcuSvF0T2NDG6qrTaiqbkFNjQkNDa3o6uxWMK5uziw4fH3c4OXlBn9/d0RG+sCLRImXK9zcnPssRJ6lo6MLJpMZra0dLKJMrR1oaTajvLyJ19bcYu52ybUAHh4uCAr24qhNz+bu/vzhND+/DosWDx3QuoSB058WczINXLNmFfbvP8THTr1z6dJlbNiwDllZWb1NLnpCxI7GmThxPLZs2Wb18ynVRR/etWtXYevWHbq+WhFeT1p6GEd3Ro+OgLOzY6I7JCZIQJSXNaGsrIkFBkVX/HzdERzsheAQLyQlBbG4cdSaCIpuURTHmvoaWjMJMtoePazi27b2ThaQISFeCA/3QViYN/+sFC0StIWLizMLHSrepTZzI9De3o7s7Fykpqbi3r370BuDLFoZ6GEnSLF3t93RQXVgV4aOZtiwZISGhuLcufP9Krij6bdbtmxHW5u2rM77i0R2XgzX7QwCRo+OtPlrUzEnCZuiwnoUFTWwIHB3c0Z4hA8iePOFl5d+xAB5l9A4jvKyZuQ+rkFxUQNHgEJDvBAT64/oaD/4+MiMLDVHdkiwrlmzGjdv3kJOTi6MhIuLC95+eyP+/OePNDLri9bYhfr6evaVexUidjQqdihi/5nPvIPNm7eyIu8P0dFRmDVrJrZu3d7v19AyIn7+coL+ZFMG3np7JB/oB/paNFmd2q1LS5v4gBkW7oPYWD/ExPhz+sco0FgOEjmJgwNRWdnSK/iaW9rh4+2G+IQAJCYGSqeWilJXdFxdtmwpHj/Ow92792BEJk2ayOfF+/cfQE9ixzhHHp2RlJSEvLy8AYmU4uISnDlzFm++uQbbtu3gri3BeFCqaNjwENy/V4H0EeF9/v66ulbk5tTgcV4tzOYujl7QSXzylDjVFz7bk8LCeowZE8m1RJTSoo2KwYnGxjael3Xq5GNuTQ8J9eLfGfnxuLo6K710w7J48SIeAWFUoUPcuHEDGza8qRGxYz0idjQKGQVSV9VAoQr806fPYt26tZ8KHmOktISnoRTW5k8ykJoW9troDkVrKGrzMLOSa278AzwweHAQli4dZqjIzesiXLRR4fSLoGjOiBHhvNHvs6qqBY9za3HjRilcnJ2QNDQIQ4cGS72PA1mwYB57zlD6ysi0t5tRUFDE7vt6SuNJGkuDaayoqEi2LT98+KjNXpNSWmQ8uG3bdrS2Gk/wSEoLPC/Lw90FI0Y+H92hw0RRYQMyMys5JUOdUBQNioz0HXAXlB6hLixKW02dFt/n76Xi5+ysGmRlV8Ns7sSQIUFITQ3TVW2TPRjI4N85c2ahtbUVFy5csumatIqnpwdWrVqJTZs2Q91IGkvXTJ48CceOnbDpa1JKizoP3nyzO8JjTSu7oC8o5UK1O6lpob1dUDU1Lci4U861JlRgO2p0JEJDvUTgvIbc3FoMG9Y/x3KK5pDgpI3ETnZ2DQ4dpHZgC1JSQzniI6ku20Fmeh0dnSJ0noCO/zU1NXxhXVKiD+NRETsadPPs6upWsraGdmoSUT0pLWlLNxZUX0NC58b1EhY7jx5Vc/s3nXRnzEwQgdMHaNDqrFkJA34dEjUpKaG8UcTnwf1KbN9+n/2GxoyNRFTUq69mhVczefJEuLm54vjxk0ovRXVcuXINM2ZM4/mKekDEjgZrda5cuWq31y8tLcORI8e4aJns0Y1gPCh0U1HRjLLSJk5VzZs/BGvXpUoEoR80N7ez2aGtxSFFfKjAmTaKuN28UYrTp/IwPCUUaWlhL60PEl7MhAlvcOrDluUAeqK2tpaNFb28vNDS0gKtIzU7GqrZIaOrjRs34P/9v4/t/l40NHTp0kWs6uvqbB9FUjNGqt8hHxwSN5Sq8vP3YCflstJGHr/wxhvRSi9Pk9y9W87OzyNHRdj9vcj1+cGDSty7V4GgQE9MmBCNgEDjzOLqb53OpEkT+LgvQufVUJEypbLOnu27l5vaanaM2xeqQcjZ0lHtgOQaumvXXqxcuRyhoTItXW9QLci1a8XsntzcbMaq1SlYvHgot0dT+/nDzCp+jtB3qKuKvHUclXqkjq6NG0dwuvH06Tzs3vUApaUy++5lTJkyCb6+viJ0rCA3NxeJiQkD9t9SAyJ2NMSIEWnIyLjrsPcjtUwztBYtWsjqXtA+NL7g/Ll8bN1yl9vEyUhw/Pjop+Y60YGNWtEpTSL0DQqUk2+OEkaB1Bm3clUKZs5MwJ075diy+S4LL4MH759i2rSp8PT0xNGjx5VeiiawWICsrGz2ddM6InY0QlhYGGpr62A2mx36vlSzs3XrNsyePRMJCX1vo9VqaPzJTQ+0tXXgzJk87Nr5gEc1kMhJTw9/6RUbFSpnZVXztHGhb3VPFB1TEkpjLVyYhGXLh3EL/ObNd1GQXwe90N/P5owZ0+Hq6iLFyH3kzp0MjBw5AlpHxI5GGDVqBG7duq3Ie5PvDg0bnTBhPM/jErQDpaIuXSzE9m33ERXliw0b05GUFPza4ll6nOp3KNUl9K3lnAwW1QD58syanYgVK4azZ8+WLXd5VpcRobE4tMufPHla6aVojubmFj4eUKGylhGxowFoR4uIiFDU78Bs7sD27TuQkjIcb7wxVrF1CNZBqYtbt0o5leHv74G33h5hlch5EjINpBlXPJVcsIrCgnrExlHDg3og0TN37hAsWZLMxdMU3aurNY6txNy5s9HZ2cFO8cJAojvp0DIidjQAFYjRYDqlIfv7Xbv2wM/PH3PnzuErJSOgtZQWOfeSOSDNqaJ0FRnR9acNmr5n8uRYXLxYaJd16g1K+Q1y6i4aViM0bX3hoqGYNj0ex4/n4sSJXE5v6jWtTLv8kiWL0NTUrOJuIm2QpYO6HXV+KoWnoHxpRkYG1AI5LdfV1WHlyhW9TruC8jQ0tGHP7kzcvVfBhapUeDzQLoqExEBUV7fwawuvpqCgDvFxAVA7ISFeWLsuDQkJAdi29R7u3CnTXRGzk5MTjzsgZ/jLl68ovRzN09XVherqarYk0SpyplI5Li4u8PLyRH29unLt16/f4Db49evfhIeH4ztPhL9AJyqqrTlw4BEmTY7BokVDbTpHafr0BJw9o3xkUe3k5NRi8BDHtJzbAqot2vjWCLQ0m7F1yz02KtTLMXPdujV48CATt2/fUXo5uuH+/QdIS0uBVhGxowFTp+xsdU6effQoC6dPn8H69etea+gk2IfKymZs/uQup5w2bEhHWJiPzd8jIsIHnZ0Wfi/h5YKzuqoFwcHaKuKkyOykybFYsDAJJ048xrmz+Zyu1irk+Lt+/VpcvXoNmZkPlV6OrigoKERsbCy0iogdlTN8+DBVf2ipaHrv3v1YtWqFIbx41FK/Q87HZ8/m87ZkaTJ3TtlzdtWMGfH8XsKLqa42IThEuwNSAwI8sHZtKoKCPbH5kwxuodfa54y6hejCiwqR1VDjqEdBX1NTo9lUlogdleedfXx8eKSFmiH/ny1btvL04PT0NKWXo3tqa0zYvDmDT1CrV6fwsE5HeLdQgSsVPwvPk51VjaQkdbSc9xcSaqmpYVixMoXTlmRZQKJaC9DoB5rnd/jwEa7TEezDvXv3kZqqzVSWiB0VEx8fh/x8bVxNkxcPmQ9GR0exAaFGL3BVf2V1+1YZDh/JxuLFyTwmwJGRhKlT43DufIHuilltQT4VJ8ervzjZGkjUrlmbyq7a5LRdX98KNRMREY5Vq5bzHL+Kikqll2OAVFY0tIiIHRVDBn6ZmY+gFegqkObNUCRqzZrVcHW1XZGs0VNa1NZMnVZNTW1Yvz6dozqOxtvbDdHRfuysLDw9goPazdXact4fSESPGRuJ+QuSsH/fQ2RnV6vy8zR4cCL76JDpqdEGFiuBxWKBydSqSYNB/Xw6dUhISDAP5NQa16/f5G6tjRvfhJ+fr9LL0TxVVS18hT1mTCSmTotXdCjfxIkxuHqlWNNFrLYmN7cGgxO1ncJ6GUFBnli/YQSyHlXj1KnHqkprjRo1EmPGjMbmzdv4BCw4huzsHCQlDYHWELGjUvz9/VTXbt4X8vLysX//Qfa6oHSc0D8ePKjEsaM5WLFyOOJUkCZxc3NGWnoYbsiQ0F5ysmswROP1Oq+CIlaLlyQjOMiLfXlaWpR31J4+fSoiIyOwY8cudHSo3xhRT+TkiNgRbNxynpOjzpZza6mpqcUnn2zGuHFjMXHiBOgZWw8PpXAxXUlTQfCb69MUmaL9MkaOjEDWoyoZI/HElHNHFIkrzYiR4Zg5KwE7tt+zmQ1BXz8zFNVcunQxOjo6cejQEakfU2hWlqenp+Y6D0XsqJTExEQ8fvwYWqe93cxXX92Opit0X8djCzo6urg+hwQO1UyozaWaTjhTp8bj/LkCGJ2SkkZERhknVRse7oNVq1M52piTU+PQ93Zzc8O6dWuRn1+AixcvOfS9hacpKytHeHgYtIS6jqICQ4LZw8NDV3loOjiRm+lbb61HYKB2XGYdTXNzO6cK0keEs3eOWolPCEBDY5tuXHf7CxVrD00KhpGgbq11b6bhbkY5rl8vcVhrOdUA0nHk7t17DnlP4eVQl7DWyhNE7KiQ4GBtFia/DjL62rVrL4ehhw7V9lA5e1Bf14qdOx5gztxEDBmi/hqQWbMSceqUsc3bSg0W2enB1dWZ68jq6kzsumzPdFJsbAxWrlyGPXv2obCwyG7vI1hPYaH23JRdlF6A8Dy0E+n1Q01t6VTHs3DhAsTFxeHkyZOq6vCwFU/WIJg7jr32+TRs8+CBLCxbPkyRtvL+dur4+bojL68WCQmOj9bZs+Xfmr8ZDUelKIeS3XFKQjUbc+cOwYXzBTh2LAfz5g15bR1HX/9mNASZXOQ3b96Ktrb2Aa5YsBX0t3B3d4OWkMiOComLI7FTCL1CxYXUqVVeXo633tpg+Lla5eVNOHQwCytXDdeM0Olh2vR4PtnpUbC+DvKeSTJYCutFTJkax51a+/c/stl+QJppzpzZbBi4det2EToqpKqqGqGh2hkdIWJHhdCIiMbGJugdyr0fPHiEQ9RGTWuVlTXh+LFcrF6TqqqOK2vx8HBB8rAQ3LlTBqORq7Ep5/Zk7LgoJMQH4MCBRwNOabm6urApaV1dHY4cOSYdVyqlqKiYHfO1wiCLwfckSqtQ8Vu37lM+HE0tfYsXL+QOJqPg4uKM+fPnobOzE8ePn9C1Yd2T6ZGqqmYcPpTN1vyentrtUqOr+U0f3+GiVRoxYC1KD1MdaGfQ6tUrOb0ykFSY3qBxJiWljVi0KKk3pdWXv3NgYACWLVuKs2fPsVeXoF4CAwMxadIEHDx4WMFVkHzpQn19/WszBBLZURkUti0tNdZVMqW16ANDdUqU1goIUN48z97U1ppY6KxclaJpofOXVvQ4nD1jnJMTmaplZWUrvQzVMWp0BMLCvHH0aE6fIzLkLbZs2RLs2bNXhI4GqK2t1VRnrYgdFYodqmUxIg8eZGLfvoNYunQRW8HrFSps3b/vEZYvH8YFrnogITGQzfVsZTanhbl1jx5pZ26dIyHLBH9/D5w9a71gmTZtCtLSUrFp0xZNO8cbDYvFopkCfenGUhnh4eG4fTsDRoXy9Js2bca0aVOxZs0qHDx4SFd+Q95ei7Fhw3ou0H7vvRpdpT3mzEnEwUNZ+Mw7fwc9QyaPlG5uanq1sLM2faP1v/vLftYFC+bh779W8dp04PLlS3maNrWWC9qiqqoKISEhmpg2L5EdlUHTZFtajG3URjUgZ86cw6VLl7F+/TqebKwHqIaBZoWdPn0GNTWOdZ91BH7+HoiJ8ecrdD0THx8vaRYrOHr0GIYMGfLSz29QUBAbBV69eo03QatOyuHQAiJ2VHbF2NWl3+LcvlJSUoqPP96MlJThWLhwPhcyaxn6GShVR1exeoWmotMsNDc3bdchvQryfcnMfKj0MlQPlezs3r0XkydPQkhI8HO/wyVLFmLnzj26/jwYIbITHKwN+wVJY6kIKszV4xX/QDCbzZzySU4eirff3oijR49rsoB7zJhRaGtrQ0bGXZsZESrJq1I0VLw7ffp07qzTG9RgRBEJW35OX/a71PI+0ANNJCfBs3btKnzyyRb+et68uRzlpHS1njsvjUB1dfVzQlatiNhREd0H0Vqll6FKHj3KYl+HRYsWcF3P6dNnuVVdC0RERCA5ORlbtmyDEcjOzsHo0aP4io8OhnoiMjKSI46C9TQ1NeHEiVPsnUPR6+vXb3KEU4AuBj27aSSKK2ksFWHrK0a9QbVM5D9EkZ133nkLUVGRUDs00HXhwnnYu3efoczRjh07gXnz5kBvDB8+HJmZcqLuK1TQHRQUyDUeInT0hcXSHfFUOyJ2VAQdDCSy83roYLlt23ZMnDgBc+fOhrOzemt5aOgpnfhbWkwwEhR9Kysr43orPREdHYnSUonsWIuTkxPmz5/LHjq/+93v2cBVCxcpQt+Mef38yJhX3YjYURF0IGhoqFd6GZqAxMPOnbtRXFzCUR412paTV1BFRQWvsT/1EE9uaqCv6zl37gLGj39DcwMDX0ZwcHea2VEBOq3vA3Q8I5NQSvuRaSilnQ8cOMg1O66u2kh9CK+n273YF2pHxI6KoKI9Iw5UHAjUFUODAqkDiMZsqOXESgf69PQ0nD9/AUaFTm5nzpzlgY56ICUlBffvP1B6GZqA9v3ly5fgwIFDuHfvfu/95JlFoyCoM1HQB/X1PSOX1I2IHUHzmEwmNiSj9BYZ9imdOqH89ZIli9gQ0ejilfxo6CpeD6mLhIQ45OeLv86roIuNlSuX8zRsso2gkQLP8vhxHndZ6sU/y+g0cBpL/ZEd6cZSCS4uLprpLlLziZXma02fPpWvLI8cOaqI9Tx1IuXmPrZp/ZWj2tLtkS45duw41q5dg48++liz4o9mANXV1Su6fiVa1PuyP8TERHMN3cmTp1/rnXPy5ClOcRUWFsJs7rDBSgWlaGhoRErKq4dwqgGJ7KgEHx8fNDY2Kr0MzUOC8dSpM3wwpeLgSZMmOnR2i5eXJwutq1evOuw9tVBflZGRwQXlWiUtTVJYL4M+XzNnTuf6LJoCb41JILUsX7hwCTNnznDIGgX70dhIkR0RO4KV+Pr6sB+FYBuqqqo5jN7a2op33nkb8fFxDnlfKr4kTxGtRjDsxa1bd5CYmKCJ3P6LSEigERF5Si9DdQQGBmDjxg2ora3jhoHW1jarv5emxtNJMiws1K5rFOxLR0enqjtie5A0lkrw9vYWsWMHbt26zb4os2bNxNixY9jVl8Ku9oA6wiiy1J/uK1ulFvqT0nBUp8+RI8e4MFVr5op0Qqd0qFoF7LN/P0ftA2+8MQ7JyUlcm0Zipz/Q53HRooUcERIEeyKRHRWZz/XlqkiwHvq9Hjp0hDujli5dwjU99rgSoVA+pdCEl0fbqBVf6QLyvpKaKimsJ/H398PGjev5M0QjH/ordAgSkTRfiXx4BMGeiNhRURcDzU4S7EdFRSUfnKura9ibh+Zt2YqhQ5PYT6S5udlmr6lHzp49j/Hjx8HT0wNaISEhgTuIhO4Zb8uWLeWJ5pcvX7GJ59C5c+cxdepktt4QBHshYkcluLu7i9hxEHSVTqInOjoaGza8ibCwsAG9Hh2kJ0+eyAWXwquhNN/x4ycxf746TPKsGc5LrbVdXcYeWEk1hevXr4Onpxc+/vgTvmCwFW1t7eyXRYX9gjbp6uriuWdqRmp2VCV22pVehmEgnw/q2KKQPNXz0If11KnTaGzse90UpWWo2LK9Xfm/n1qcdl8F1TSlpaXydHQaGqpmtJjCsvU+QE7gI0ak4/DhI6isrII9uHHjFt5+ewPu3r3rMIdqwXZQI4i7uwfPL1Qr6pZiBsLNzQ3t7RLZcTRUM7B7915cu3adw/OzZ8/kv0VfoMLn69dv2G2NeoSE5pQpk1XjeP0yqIPMqCks8hai2hxqnvj44012EzpER0cHXzCQS7WgPTo6OlTfkSViRyWQy6yYaykHTVKn1FZRUTE2bnwT48aNscqfhworyciQfEME66F9/fTpM9yqr1aoTZ68r4yWwqL9nmpoqHOOojkXLlx0SCfatWs3+HMnaDM97axysTPIYjF20JCGmFFenpw8lTRG2rt3L+bPn89dWYKy0Mnt9u3ulvUxY8Zg2LBhLy2e3LZtG5YsWQIvLy+Hr1MPHDlyBMnJyVwErDYuXbqE0NBQDBkyBEaBJtWfOnUKaWlpSE9Pd3jR8OnTp3lfiI+Pd+j7CgPj+PHjfKwMCgqCI6F6utjYWNTV1b3Ww8vwNTs9rsX0CxMEQRAEQXvn8deJHcNHdugqvqSkBL6+vtL6KDzF0aNH8ZWvfIWvdqdMmYIPPvgAkZHaH2ipNubMmcNRtIKCAp4Rpxbee+89HDt2jKO+RuBXv/oV/sf/+B/cLPHb3/4Wq1atUnpJgvBKSL6Q0ImKioKT06urcgwvdgThdfyf//N/8A//8A/cbUU+Oq/7UAl9g4wGyQZg/fr1+Oijj6CmAt25c+dyqlLvvPHGG7h+/TrWrl2Ljz/+uM9F+oKgdkTsCIKVEUAa7jlx4kSll6JLvv/97+Of//mfcevWLYwcOVLp5XAL9IgRI/hvTkJA7/ziF7/A1KlTZf8WdIuIHUEQVAGJiy9+8YucOlSaz3zmM9izZw8XQAqCoH1E7AiCIDwDiZyMjAyOdgiCoH1E7AiCIAiCoGuk0lIQBEEQBF0jYkcQBEEQBF0jYkcQBEEQBF0jYkcQBEEQBF0jYkfQNWfOnMHy5cvZYZMcsnft2vXU41Sf/73vfY+dkT09PTFv3jxkZWU99Zyamhq88847PDuN5qh97nOfQ1NT01PP+f3vf8/zfGg+zOXLl2F0yDOHft9PbsOHD+99vLW1FV/+8pcRHBwMHx8fNrMrLy9/6jWo9ZvmZtFssn379tl9zb/+9a95LhPNpyO/mStXrvQ+9vDhQ+7MiomJwQ9+8AOoHXJAJr8i2mdpmzx5Mg4ePPiUUeasWbP4Mfrb0GyhZ5H9XtAV1I0lCHrlwIEDlu985zuWHTt2UNehZefOnU89/pOf/MTi7+9v2bVrl+X27duWFStWWBITEy0mk6n3OYsWLbKMGjXKcunSJcvZs2ctSUlJlrfeeqv38fz8fL7vwoULlq1bt1pSUlIsRuf999+3pKWlWUpLS3u3ysrK3sf/23/7b5bY2FjL8ePHLdeuXbNMmjTJMmXKlN7HW1tbLTExMZajR49ajhw5wv9ua2uz23o/+eQTi5ubm+WDDz6w3Lt3z/KFL3zBEhAQYCkvL+fH582bZ/nNb37Da33jjTcs58+ft6iZPXv2WPbv32959OiR5eHDh5Zvf/vbFldXV8vdu3f58V/+8peWH//4x7zR56K2tva515D9XtATInYEw/Cs2Onq6rJERERY/uVf/qX3vrq6Oou7u7tl06ZN/PX9+/f5+65evdr7nIMHD1oGDRpkKS4u5q8zMjL4BNjU1GTJzc21JCQkWIwOiR06Ub4I+h3TiZdOkD08ePCAf88XL17kr+vr6y3x8fEskGij32lDQ4Pd1jthwgTLl7/85d6vOzs7LVFRUSwGiHHjxlkuX75saW9vZ0FMQkJrBAYGWv7v//2/T9138uTJF4od2e8FvSFpLMGwPH78mId8UuqqB5qcSymMixcv8td0SyH8J0cG0PNpPlZP2D49PZ1TBvS9aWlpmkhzOAJKB1L6cPDgwZwOoUGfBM1gMpvNT/3eKcUVFxfX+3un1AkN4qT0Ir3Gl770JR7Waw9o5hmt6cn10N+Xvu5ZD42zoK+9vLz4sYULF0IrdHZ24pNPPuG5bpTOsgbZ7wW9oZ4Rw4LgYEjoEOHh4U/dT1/3PEa3YWFhTz1Ok7mDgoJ6n0P84Q9/wE9/+lM+GVLtj9EhwfjHP/6R621KS0t5mvb06dN55hT93mjQJJ1MX/Z7J95//3387d/+LZ9g7SV0iKqqKhYEL9oPaBo7sWTJElRWVrKzcmhoKLQAOUCTuKH6KKqL2rlzJ1JTU636XtnvBb0hYkcQbAQV2wrdLF68uPffdPVP4ocKWbds2dKnkyJFDdSCu7u7ZoQOQUKTBqvW19fz5PZ3330Xp0+ftlrwWIvs94IWkDSWYFgiIiL49tkuIPq65zG6raioeOrxjo4O7lTpeY7weiiKQ51V2dnZ/Huj1NGzHUBP/t4dSUhICJydnV+5H2gRip4lJSVh3Lhx+PGPf4xRo0bhV7/6lVXfK/u9oDdE7AiGJTExkQ/cx48f772P0hRUk9BT20C3dFKmmo4eTpw4ga6uLo5WCNZBLcs5OTlcg0MnX1dX16d+79TaTTU91taU2FoU0JqeXA/9felrJdZjL+hnamtrs+q5st8LekPSWILuT7IUTXiyKJlC+1R7QAWxVBNChZVDhw5l8fPd736XC2JXrVrFz09JScGiRYvwhS98Ab/73e+4sPYrX/kKNm7cyM8TXszXv/519jei1FVJSQnX31D05K233uLUFHm2fO1rX+O/AxUjf/WrX+UT7KRJkxRZL62F0jxUkDthwgT867/+Kxf0UpG0FvnWt77FqUTaxxsbG/Hxxx/j1KlTOHz4MD9OdTe09Xw2qL6H6qLo+fQ3kf1e0B1Kt4MJgj3paa19dnv33Xd728+/+93vWsLDw7nlfO7cuexL8iTV1dXsL+Lj42Px8/OzvPfee5bGxkaFfiJtsGHDBktkZCR710RHR/PX2dnZvY+Tj9Hf/M3fcDu0l5eXZfXq1ezFoyT/9m//ZomLi+M1Uys6+ctolb/+67/m1n36WUJDQ3m/Jr+iJ60BXvS5+PDDD3ufI/u9oCcG0f+UFlyCIAiCIAj2Qmp2BEEQBEHQNSJ2BEEQBEHQNSJ2BEEQBEHQNSJ2BEEQBEHQNSJ2BEEQBEHQNSJ2BEEQBEHQNSJ2BEEQBEHQNSJ2BEEQBEHQNSJ2BEEQBEHQNSJ2BEHQFX/84x8xaNAg3mj22atISEjgOVi25J//+Z9739/Wry0IQv8QsSMIgt357Gc/2ysAaOI5DV39x3/8R7S2tlr1/QsXLuRBolevXrXq+TRctLS0FP/zf/5PKDEEld47JibG4e8tCMKLkanngiA4BJqi/eGHH/IE7evXr/OUcRI//+t//a9Xfl9BQQEuXLjAU7c/+OADjB8//rXvRa8bEREBJfDx8eGNxJkgCOpAIjuCIDgEd3d3FiCxsbFYtWoV5s2bh6NHj772+0ggLVu2DF/60pewadMmmEymfr1/RUUFli9fDk9PT44sffTRR889p66uDp///OcRGhrK0aE5c+bg9u3bTz3nBz/4AcLCwuDr68vP/eY3v4nRo0f3a02CIDgGETuCIDicu3fvcrTGzc3tlc+zWCwsdj7zmc9g+PDhSEpKwrZt2/qdSissLMTJkyf5NX7zm9+wAHqSN998k+87ePAgR5/Gjh2LuXPnoqamhh8ngfTDH/6Qo1H0eFxcHH7729/2az2CIDgOSWMJguAQ9u3bx+mdjo4OtLW1wcnJCf/+7//+yu85duwYWlpauGaHINHzhz/8Af/lv/yXPr33o0ePWMBcuXKlNw1Gr5OSktL7nHPnzvHjJHYoCkX87Gc/w65du1gc/df/+l/xb//2b/jc5z6H9957jx//3ve+hyNHjqCpqanPvw9BEByHRHYEQXAIs2fPxq1bt3D58mWu1yHBsHbt2ld+D9XobNiwAS4u3ddlb731Fs6fP4+cnJw+vfeDBw/4NcaNG9d7H0WKAgICer+mdBWJluDg4N66G9oeP37c+34PHz7EhAkTnnrtZ78WBEF9SGRHEASH4O3tzWmoHhEzatQojq5QpORFUOpo586dXND8ZKqos7OTv5/SSbaEhE5kZCROnTr13GNPiiJBELSHRHYEQXA4lML69re/jX/6p396acEx1cdQ+zZFXCgi1LP9/Oc/Zy8dEj3WQlEcSp9RnU0PFKWhguQeqD6nrKyMI0Akyp7cQkJC+DnDhg17rv3d2nZ4QRCUQ8SOIAiKQMXA1J7961//mr+mKA6Jkh4o6rNu3Tqkp6c/tVEkqKqqCocOHbL6vUikUOv7F7/4RU6jkeihTirqzOqBusMmT57MnWJUh5OXl8dF1N/5zndw7do1fs5Xv/pVXtef/vQnZGVlcWfWnTt3uNVdEAT1ImJHEARFoAgKeef89Kc/RXNzM+rr6znaQpAYoYjOi2p6/P39uUOKREdfoK6uqKgozJw5E2vWrOGCY2oh74EEy4EDBzBjxgyuJ0pOTsbGjRuRn5+P8PBwfs4777yDb33rW2wcSJEgquehLi8PD48B/z4EQbAfgyzU2ykIgqATKMVFYyKeTFHZk/nz57N/0J///OfnRlHQOl43skIQBPsjkR1BEHQHRYmok+ob3/iGTV+X2uB/8Ytf4N69e8jMzMT777/P7fHUXdbDj370I35vcn4WBEEdSGRHEARd0djYiPLy8t4uqp7iYltAxdTkwnzz5k2e60W1QFRkTWmxJ7vIekwIyYmZ0m6CICiLiB1BEARBEHSNpLEEQRAEQdA1InYEQRAEQdA1InYEQRAEQdA1InYEQRAEQdA1InYEQRAEQdA1InYEQRAEQdA1InYEQRAEQdA1InYEQRAEQYCe+f8WXAacjJtibgAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_map(ps.W_HI,ps.wproj,have_cbar=False)" ] }, { "cell_type": "code", "execution_count": 5, "id": "0c9b76d4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "total area of the survey is 7349.0 sqdeg\n", "redshift range is 0.5 - 1.5\n" ] } ], "source": [ "print(\"total area of the survey is %.1f sqdeg\" % (ps.W_HI[:,:,0].sum() * ps.pixel_area))\n", "print(\"redshift range is %.1f - %.1f\" % (ps.z_ch[-1], ps.z_ch[0]))" ] }, { "cell_type": "markdown", "id": "9187bad0", "metadata": {}, "source": [ "## Box-dimensions and k-modes\n", "\n", "You can then retrieve the lightcone dimensions very easily:" ] }, { "cell_type": "code", "execution_count": 6, "id": "373febd2", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([5745.6529501 , 7940.18852805, 3781.26218245])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ps.get_enclosing_box()\n", "ps.box_len" ] }, { "cell_type": "markdown", "id": "165e73ef", "metadata": {}, "source": [ "this box dimension is a bit arbitrary. Basically you just find a box that encloses all the pixels in the map cube. There are some choices you need to make, such as the buffer at each side (`ps.box_buffkick`), the resolution relative to the map (`ps.downres_factor_radial`, `ps.downres_factor_transverse`) etc. If you have the exact box dimensions from Isa you should just manually put them in. For example:\n", "\n", "```python\n", "ps._box_len = np.array([200,600,300]) # just for illustration\n", "ps._box_ndim = np.array([20,60,150]) # just for illustration\n", "\n", "# if you manually put in you also need to do this line\n", "ps.propagate_field_k_to_model()\n", "```" ] }, { "cell_type": "markdown", "id": "114817e1", "metadata": {}, "source": [ "Note that, the k-sampling also depends on your choice of grid resolution. Two factors come into play here: First, of course the grid resolution sets the Nyquist frequency of the k-grid; Second, the compensation effect of gridding is determined by the Nyquist frequency as well (see cookbook about gridding compensation). To ensure completeness in the sampling, it is typical to choose a resolution lower than the map pixel and channel resolution. This is controlled by these two parameters:" ] }, { "cell_type": "code", "execution_count": 7, "id": "892aefae", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([66.04198793, 66.72427334, 5.00829428]),\n", " array([5745.6529501 , 7940.18852805, 3781.26218245]))" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# downscale the resolution along line-of-sight\n", "ps.downres_factor_radial = 1.5\n", "# downscale the resolution on the transverse plane\n", "ps.downres_factor_transverse = 1.2\n", "ps.get_enclosing_box()\n", "# note how anisotropic the box is\n", "ps.box_resol,ps.box_len" ] }, { "cell_type": "markdown", "id": "6eb3e169", "metadata": {}, "source": [ "If you want a specific resolution in Mpc, you can check the pixel and channel resolution first and calculate these scaling factors:" ] }, { "cell_type": "code", "execution_count": 8, "id": "591fcaa1", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(np.float64(55.55517205880652), np.float64(3.337730934154746))" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ps.pix_resol_in_mpc, ps.los_resol_in_mpc\n", "# set scaling factor \n", "# ps.downres_factor_transverse = target_res_pix / ps.pix_resol_in_mpc\n", "# ps.downres_factor_radial = target_res_los / ps.los_resol_in_mpc" ] }, { "cell_type": "markdown", "id": "4c9023ec", "metadata": {}, "source": [ "Setting the box dimensions gives you the corresponding k-modes:" ] }, { "cell_type": "code", "execution_count": 9, "id": "64ca2fb6", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "((87, 119, 378), (87, 119, 378), array([0.04702285, 0.04668755, 0.62644713]))" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# |k| and \\mu=k_para/|k|, with the 3D shape of the box\n", "ps.kmode.shape,ps.mumode.shape, ps.k_nyquist" ] }, { "cell_type": "markdown", "id": "9471200d", "metadata": {}, "source": [ "Given the k and mu mode, you can then propagate the model PS to the k-grids. Here we will just use a simple linear matter power spectrum with Kaiser RSD. You can take the k-values and mu-values to be input for your own model, and add in observational effects such as the beam etc (or use `meer21cm`, see other tutorials and here we focus on just highlighting the k-sampling) " ] }, { "cell_type": "code", "execution_count": 10, "id": "2895e8b9", "metadata": {}, "outputs": [], "source": [ "power_model_3d = ps.auto_power_matter_model" ] }, { "cell_type": "markdown", "id": "bdaed0f0", "metadata": {}, "source": [ "Suppose you bin the monopole to 1D:" ] }, { "cell_type": "code", "execution_count": 11, "id": "d92d0b30", "metadata": {}, "outputs": [], "source": [ "ps.k1dbins = np.linspace(0, 0.5, 21)\n", "power_model_1d_1, keff, nmodes = ps.get_1d_power(power_model_3d)" ] }, { "cell_type": "code", "execution_count": 12, "id": "31b4739e", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAG1CAYAAADwRl5QAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPRpJREFUeJzt3Qd0VVX6/vEnPSQkoSeEBALSDCWEELpjGRxAZRSUQRFpIzbUcXQsWH8zf3UcHRtj7BQRUZSmqNiwIZ2EJh0hJNTQUiCk819nYyII6AWSnHvP/X7Wust7Ty7JGwHzuM9+3+1z9OjRowIAAHAgX7sLAAAAqCoEHQAA4FgEHQAA4FgEHQAA4FgEHQAA4FgEHQAA4FgEHQAA4FgEHQAA4Fj+8nJlZWXauXOnwsLC5OPjY3c5AADABda847y8PEVHR8vX9/TrNl4fdKyQExsba3cZAADgLGRmZiomJua0H/f6oGOt5JT/iwoPD7e7HAAA4ILc3FyzUFH+c/x0vD7olN+uskIOQQcAAM/ye9tO2IwMAAAci6ADAAAci6ADAAAci6ADAAAcy2uDTkpKiuLj45WcnGx3KQAAoIr4HLUm7nh5e1pERIRycnLougIAwGE/v712RQcAADgfQQcAADgWQQcAADgWQQcAADgWQQcAADgWQaeKeHkzGwAAbsHrD/WsKo98+KP25RVpRI84dW5a53cPHQMAAJWPoFMFDhWWaFrqdhUUl+mzNbvVJjpcI3o0Vb+Ehgry97O7PAAAvAYDA6toYODGPXmaMD9dM5cfCzyWejUDNbhLEw3p2lgNwoIr7WsBAOBtcl38+U3QqeLJyAcPF+m9pZmatDBdu3IKzLUAPx/1ax9tVnnaxURU+tcEAMDpcgk67nUERHFpmT5fs9us8qRuO1hxvVOT2ibw9G4TKX8/9oYDAOAKgo4bn3W1MjNbE+Zv1Serd6m49Ni//uiIYA3tHqdrk2NVKySwWuoAAMBTEXQ84FDPrNwCTV60Te8sztD+w0XmWnCArwZ0jNGI7nFqERlWrfUAAOApCDoedHp5QXGpPlq509zWWrcrt+L6BS3qaWSPprqwZX35+tKeDgBAOYKOBwWdctZvxeKtBzT+h636ct0elf/ONKsXquE94nR1xxiFBjERAACAXILOb0tJSTGP0tJSbdy40S2CzvEyD+TrrQXpmro0U3mFJeZaWLC/BnWK1bDucYqtE2J3iQAA2Iag44ErOqcbPjg9dbsmLkjX1n2HzTXrLtal8ZGmW6sLU5cBAF4ol6DjjKBTrqzsqL7buFfj52/VvE37Kq6f3zBcI3vEqV9CtIIDmLoMAPAOuQQdZwWd422ypi4vSNeMtF+mLtcNDdT1XRprSNcmahDO1GUAgLMRdBwcdMpl5/88dXlBunYeN3X58nYNzW2thNhadpcIAECVIOh4QdApV2KmLu8xQwiXHTd1OclMXY5TnzZRTF0GADgKQceLgs7xVm/PMYFn9qqdFVOXG0YE64ZuTXRdcmPVDmXqMgDA8xF0vDTonDB1eXGGpizepn2Hfpm63D8xxqzytGTqMgDAgxF0vDzoHD91+eNVu8wQwrW/mrpsBZ6LWjZg6jIAwOMQdFzk9KBTzvptXrL1gDlm4ou1u1X28+9603qhGtatia7pFKuaTF0GAHgIgo6LvCXo/Hrq8qSF6aZjK6/g56nLQf76S3KshnWLU+O6TF0GALg3go6LvDHolDtcWGJm8VirPFt+nrpsDVnudb41dTlO3ZrVZeoyAMAtEXRc5M1B54Spy5v2msDz/ca9FddbR4WZ09P/3IGpywAA90LQcRFB50Sbs/JM4JmRtkNHikvNtTqhgRrcubFpUY9k6jIAwA0QdFxE0Dm1nPxivbc0Q5MWbtOO7CPmmr+vjy5vf2zqcgemLgMAbETQcRFB5/enLn+x9tjU5aXpv0xdTmxcy9zW6tM2SgFMXQYAVDOCjosIOuc2dTkq/NjUZevWFlOXAQDVhaDjIoLOmcvKK9A7izL0znFTl4P8fTWgYyMN795UraKYugwAqFoEHRcRdM5eYUmpZq/cZVZ51uz8Zepyj+Z1zW2ti1sxdRkAUDUIOi4i6Jw764+QtX/HCjyfr/ll6nJc3RAN6x6ngUxdBgBUMoKOiwg6lT91+e1F2/Tekgzl/jx12Qo5AzvFaHj3ODWpG2p3iQAAByDouIigU4VTl5fvMKs8W/b+MnX5j60jNdKaunweU5cBAGePoOMigk7VT13+/uepy9/9auqytcJzVWIjpi4DAM4YQed3pKSkmEdpaak2btxI0KkGm7MOaeKCrZqe+svU5dohARrcpbFu6BqnqAimLgMAXEPQcRErOvZMXZ66LENvLThx6vJl7aypy3FKbFzb7hIBAG6OoOMigo69U5e/NFOX07Uk/UDFdet4CSvwWMGHqcsAgFMh6LiIoOMeftyRo/Hzt+rjlbtUVFp2wtTl6zo3NgeLAgBQjqDjIoKOe9mbV2gmLk9elKF9hworpi5f1aGRRvSMU+sofo8AACLouIqg475Tl63VnQkLturHHb9MXe5+Xl1zevolrRvIj6nLAOC1cgk6riHouDfrj+eybQc1/ocTpy43rnNs6vJfOsUoLDjA7jIBANWMoOMigo7n2H4wX28v3KZ3fzV1+ZqkY1OX4+oxdRkAvEUuQcc1BB3Pk19Uohlpx6Yu/3Tc1OVLWjXQyJ5Nze0tpi4DgLMRdFxE0PHsqcvzNu8zgefbDb9MXW4VGabhPeLUn6nLAOBYBB0XEXSc4ae9hzRxfrqmp21XftGxqcu1rKnLnRubFvWGETXsLhEAUIkIOi4i6DhLzpFivb80UxMXpFdMXba6s/q2jTK3tToydRkAHIGg4yKCjnOnLn+1bo/GW1OXt/4ydTkhtpY5Pb1v24YK9GfqMgB4KoKOiwg63jF12Vrh+WjFzoqpyw3CgjT056nLdWsG2V0iAOAMEXRcRNDxrqnLUxZnaPLibea5xVrV6dc+Wtd2jlWnJrXp1gIAD0HQcRFBx/sUlZTpk9U7Nf6HdK3ekVNxvVn9UF2bHKsBHWNUj1UeAHBrBB0XEXS8l/VHPy3joN5bkqmPV+3SkeJj3Vr+vj66ND5Sg5JjdUGL+hw1AQBuiKDjIoIOLHkFxZq9cpemLsvUyszsiuvREcG6plOsBibFKLZOiK01AgB+QdBxEUEHv7ZuV66mLs3UzOU7TLu6xdq607N5PV2b3Fi94hsoyJ9BhABgJ4KOiwg6OJ2C4lJzkKgVehb8tL/iep3QQA1IbGRubbWIDLO1RgDwVrkEHdcQdOCKbfsP6/1lmfpg2XZl/dyxZenYuJZZ5bm8fUOFBvnbWiMAeJNcgo5rCDo400GE323cq/eWZurr9VkqLTv21yc00E9/7hCtQcmNlRATQZs6AFQxgo6LCDo4W1m5BZqWtt3c2tq2P7/ieuuoMHNbyzpUtFZIoK01AoBTEXRcRNBBZZyivnjrAU1dmqE5P+5WYUlZxTDCPm2izGyers3qypc2dQCoNAQdFxF0UJly8ov14codendJpuneKhdbp4YGdYrVNUmxiooItrVGAHACgo6LCDqoCtZfqx935Oq9pRnmjK28whJz3VrUubhVA3Nr6+LWDRTgx8GiAHA2CDouIuigquUXlejT1VabeoaWph+suF4/LEhXd4wxoadpvVBbawQAT0PQcRFBB9Vpc9YhfbAsU9NSt2v/4aKK612a1jEHi/Zt21DBAQwjBIDfQ9BxEUEHdh0s+vX6PaZN3WpXL/9bGBbsb7q1rFWeNtERdpcJAG6LoOMigg7stiP7iKYt224GElrPy7VrFGECjzWfJzw4wNYaAcDdEHRcRNCBO7Wpz/9pn1nl+WLNbhWXHvurGRzgq8vaNTQTmJPjajOMEABE0HEZQQfuaP+hQnOoqDWMcFPWoYrrzeqFmlWeAR1jzGZmAPBWuQSd35aSkmIepaWl2rhxI0EHbsn665mWkW06tmav3KUjxaXmur+vj3qdH6lBnWP1hxb15ccwQgBeJpeg4xpWdOAp8gqK9fGqXebW1srM7IrrDSOCNbBTrAYmxSi2ToitNQJAdSHouIigA0+0fneuua1l3d7Kzi8216ytOz2b1zO3ti6Nj1SQP23qAJyLoOMigg48WUFxqb5Yu8fc2pq/eX/F9dohAWYfjxV6WkaG2VojAFQFgo6LCDpwioz9+fogNdO0qe/JLay43rFxLdOxdXn7hgoN8re1RgCoLAQdFxF04DQlpWX6ftNevbckU3PXZ6m07Nhf8dBAP/VLiDb7eazwQ5s6AE9G0HERQQdOlpVXoOmpVpt6htL351dct87WGpDYSFclNmIDMwCPRNBxEUEH3sD6a7546wG9vzRTc37cXdGmbunarI7Zz2MNJazJrS0AHoKg4yKCDrzNocISffbjbs1I266FW/ZXnLNlTWDu0ybKhJ4ezesxmweAWyPouIigA29mna01a/kOTU/dri37DldcjwwPMre1ru4YQ9cWALdE0HERQQc4dmtrRWa2ZqTt0EcrdyrnyLHZPJa2jcJN4PlzQrTq1uTYCQDugaDjIoIOcKLCklJ9sz5L09N2mH+W/Ny1ZR07cVGrBrq6YyNdcn4DBhICsBVBx0UEHeC3DxedvXKnZizfoVXbcyquR9QI0BXtG+rqpBglxtKqDqD6EXRcRNABXLNpT55Z5bH29OzOLTipVb1/x0aKqU2rOoDqQdBxEUEHODPWAMKFP+03XVunalW39vP0pVUdQBUj6LiIoAOcPVrVAdiFoOMigg5QPa3q13SMUQta1QFUEoKOiwg6QPW1qrdrFKEBHRvRqg7gnBF0XETQAaoOreoAqgpBx0UEHcD+VvV+CQ3Nfh5a1QG4iqDjIoIO4D6t6s2sVvWOx05Vp1UdwG8h6LiIoAPY26q+4Kd9Zj/PZ79qVe/WrK4JPbSqAzgVgo6LCDqAe7WqW11bVqt6uRoBfurT1mpVb6Tu59GqDuAYgo6LCDqA+9l+MF8frth5Uqt6VHiw/u/PbUzwAeDdcgk6riHoAJ7Xqv63P7YwD19WdwCvlUvQcQ1BB/CcVvX/zNmg8fO3mtfW5OVn/5KgUPbvAF4p18Wf377VWhUAnCVr1s6j/eL19DXtFejnq8/W7NbVryxQ5oF8u0sD4MYIOgA8yl86xerdm7qqXs0grd+dpz+/9IM5ZBQAToWgA8DjJDWprdl39DBHShzML9YN4xbr7YXpZk8PAByPoAPAIzWMqKEPbummKztEm6MlHvlwjR6c+aOKSsrsLg2AGyHoAPBYwQF+emFQBz3Qt7WskyPeXZKhIW8u1r5DhXaXBsBNEHQAeDTrbKxbLjxP44clKyzIX0vSD+jKl+Zrzc5fztMC4L0IOgAc4eLWDTRzdA81rReqHdlHTEfWJ6t22V0WAJsRdAA4RvMGNTXrth76Q8v6Kigu0+gpaXr2iw0qK2OTMuCtCDoAHCUiJEAThidr1AVNzev/fb1ZN09ONWdpAfA+BB0AjmMd/PnQ5fF6dmCCAv199eXaPRrw8nxt2//LuVkAvANBB4BjXZ0Uo/dv7qYGYUHauOeQrkyZr/mb99ldFoBqRNAB4GgdYmtp9h09lRBbS9n5xRo6fokmzN/KcEHASxB0ADheZHiwpt7UVQMSG6m07Kj+OXutHpi+2hwUCsDZCDoAvGa4oHXa+cOXny9fH2nqskwNfmOx9uYxXBBwMoIOAK8aLnjjBc00YURnhQX7K3XbQXMo6OrtDBcEnIqgA8DrXNiyvj4c3UPN6odqV06Brnl1gT5cscPusgBUAYIOAK/UrH5NzRrdQxe3qq/CkjL97b0V+s9n680eHgDOQdAB4LXCgwP05rBkc1aW5ZVvf9KoScuUV1Bsd2kAKglBB4C8fbigdfr5i9d2UJC/r75en6X+Ly/Q1n0MFwScgKADAJKu7NBIH9zSTVHhwdqcdUhXvvSDvt+41+6yAJwjgg4A/Kx9TC19dHsPdWxcS7kFJRo+YYnenLeF4YKAByPoAMBxGoQH692bumpgUoysfcmPf7JO//hglQqKGS4IeCKCDgD8SpC/n56+pr0evSLe7OGZnrZd176+SFm5BXaXBsBbg05+fr6aNGmif/zjH3aXAsAhwwVH9myqt0Z0VkSNAK3IzFa/l37Qysxsu0sD4I1B54knnlDXrl3tLgOAw/RsUc8MF2zRoKb25BZq4GsLNXP5drvLAuBNQWfTpk1av369+vbta3cpABworl6oZtzWXb3Ob6CikjL9fepKPfnpOoYLAh7A9qDz/fffq1+/foqOjjZLxbNmzTrpPSkpKYqLi1NwcLC6dOmiJUuWnPBx63bVv//972qsGoC3CQsO0Os3dNLtFzc3r1//fovpyjpwuMju0gC4c9A5fPiwEhISTJg5lalTp+ruu+/WY489prS0NPPe3r17Kysry3z8ww8/VMuWLc0DAKqSr6+P/tG7lV4anKgaAX6at2mf+v2PfTuAO/M56kYDIqwVnZkzZ+qqq66quGat4CQnJ+ull14yr8vKyhQbG6s77rhDDzzwgMaMGaPJkyfLz89Phw4dUnFxse655x49+uijp/wahYWF5lEuNzfXfL6cnByFh4dXw3cJwAk27M7TLZNTzQTlQD9f/fPKNro2Odb8dwxA1bN+fkdERPzuz2/bV3R+S1FRkVJTU9WrV6+Ka76+vub1woULzWvrllVmZqbS09P13//+V6NGjTptyCl/v/UvpvxhhRwAOFOtosL04e099Kf4SBWVlmnMjNW6fzrzdgB349ZBZ9++fSotLVVkZOQJ163Xu3fvPqvPaa0AWemv/GGFJAA420NBXx2SpPv6tJKvj/T+su265tUFyjyQb3dpAH7mLwcZPnz4774nKCjIPACgsvbt3HZRcyXE1NId7y7XjztyzbydFwZ10EWtGthdHuD13HpFp169embvzZ49e064br2OioqyrS4A+LUezevp4zt6KiG2lrLzizVi4lKNnbtJZbSgA7Zy66ATGBiopKQkzZ07t+KatRnZet2tWzdbawOAX4uuVUPv39xV13dpLKvN47kvN+rGScuUk19sd2mA17I96FidUitWrDAPy9atW83zjIwM89pqLX/jjTf01ltvad26dbr11ltNS/qIESNsrhwATn1O1hP92+mZa9oryN9XX6/PMrey1uzMsbs0wCvZ3l7+7bff6uKLLz7p+rBhwzRx4kTz3Gotf+aZZ8wG5A4dOmjs2LGm7bw629MA4Ez9uCNHt76TqswDR0zoebJ/O12dFGN3WYAjuPrz2/agYzeCDoCqlJ1fpLumrtC3G/aa10O6NtYjV8SblR8AXj5HBwA8Xa2QQI0flqy7erWQNUtw8qIMDXptkXblHLG7NMAreG3QsY6ciI+PN1OXAaCqW9Dv6tVS44cnK6JGgFZkZuuKsT9oweZ9dpcGOB63rrh1BaAaZezPN0dHrN2Va4YM3tentW7+QzOOjgDOELeuAMANNa4bohm3ddfVHWNkjdh5as56E3zyCmhBB6oCQQcAqllwgJ/+O7C9nujf1hwI+vmaPbrypfnauCfP7tIA7751ZQ3r++677zRv3jxt27ZN+fn5ql+/vhITE81Bm554QCa3rgDYydqvc9vkVO3MKVBIoJ/+c3V79UuItrsswLtuXR05ckSPP/64CTKXXXaZ5syZo+zsbHM8w+bNm/XYY4+padOm5mOLFi2qzO8DABytQ2wtzb6jp3o0r6v8olJzXta/Zq9VcWmZ3aUB3rOiYwUc68gF69DMSy+9VAEBASe9x1rhmTJlil577TU99NBDGjVqlDwBKzoA3EFJaZme/XKjXvn2J/M6Oa62UgZ3VIPwYLtLA5w/MNA6euH888936QsXFxeb4xvOO+88eQKCDgB38vma3frH+yuVV1ii+mFBevn6jkqOq2N3WYDbYTKyiwg6ANzNlr2HTCfWxj2H5O/rowcvO18jesTRgg5UdXu5lYmsQzdLSkrM66KiIk2dOlWTJk3Svn2eNfiKgYEA3FWz+jU1a3QP/TkhWiVlR/Wvj9fqb++tUGFJqd2lAR7H5RWdDRs2qHfv3srMzFSzZs30xRdfaODAgVq/fr0JQCEhIVqwYIFatGghT8KKDgB3Zf23deKCdD3xyToTeHq3iTT7dvz9mAwC5Fb2is7999+vhIQErVixQldccYUuv/xyxcTE6ODBgzpw4IDZrPyvf/2rsuoHAK9n3aoa0aOpJoxIrpi3c9+0VSqzJg0CqNwVnQYNGphVnA4dOujw4cMKCwvT999/r549e5qPW6s51113nem+8iSs6ADwBF+u3WP27ZSWHdUNXZvoX1e2Yc8OvFpuZa/oHDp0SHXqHNv5Hxoaah4NGzY8oQV9z54951o3AOAULo2P1HN/STAnoL+9aJue+XyD3SUBHsHloBMdHW3axss9/fTTZpWn3N69e1W7du3KrxAAYFzZoZEev6qtef7ytz/p5W83210S4JygYx3xYG08Lnfrrbea21flrNtaHTt2rPwKAQAVru/SRGP6tjbPn/5sg95emG53SYBbq7Q5OlbbeXBw8Am3szwBe3QAeKJnv9ig/319bEXHuqU1oGOM3SUBbvnz27+yvqB11hUAoHrcfWlL5RWUmPbze6etUkigv/q0jbK7LMDtnPEwhjvvvFNjx4496fpLL72ku+66q7LqAgD8Bqvj6tEr4nVNUozpxLrz3eX6fuNeu8sCPD/oTJ8+XT169Djpevfu3TVt2jR5CiYjA/B0vr4+empAO13WLkpFpWW66e1lWpZ+wO6yAM8OOvv37zf3xH7Nuj/mScdAjB49WmvXrtXSpUvtLgUAzpo1JfmFQYm6sGV9FRSXacSEpfpxR47dZQGeG3SaN2+uzz777KTrc+bMMUdDAACqV6C/r14dkqTOcXXMqedDxy/R5qw8u8sC3MIZb0a+++67dfvtt5u5OZdccom5NnfuXD377LN64YUXqqJGAMDvqBHop3HDO2nwG4u1ekeOhry5RB/c0k2xdULsLg3wvPbyV155RU888YR27txpXsfFxen//u//NHToUHka2ssBOMmBw0Ua9NpCbco6pMZ1QkzYiQwPtrsswLaf3+c0R8da1alRo4Zq1qwpT0XQAeA0e3ILNPDVhco4kK8WDWpq6s3dVCc00O6yAPc+6+rXsrKyzKTk5cuXm8ADAHAP1grOOzd2UWR4kFnZGTZ+ifIKiu0uC7DFGQedvLw83XDDDebsqwsvvNA8rOdDhgwxqQoAYD9rb44VdqyVHGvPzl8nLtORolK7ywLcP+jceOONWrx4sT755BNlZ2ebx8cff6xly5bp5ptvrpoqAQBnrHmDME0a2VlhQf5akn5At0xOVVFJmd1lAdXqjPfohIaG6vPPP1fPnj1PuD5v3jz16dNHhw8flidhjw4Ap7OGCN4wbomOFJea4YJjr00083cAT1Zle3Tq1q17yoGB1rXatWufeaUAgCrVKa6OXrshSYF+vvp09W49MGO1ysoq5TxnwO2dcdB5+OGHzSyd3bt3V1yznt9777165JFHKrs+AEAl+EPL+hp7XaL8fH00LXW7/vXxWp1D0y3g3FtXiYmJ2rx5swoLC9W4cWNzLSMjQ0FBQWrRosUJ701LS5O749YVAG8yPXW77vlgpXl+5yXNdfefWtldElClP7/PeDLyVVddJSewDvW0HqWldCEA8B5XJ8XocFGJHv1wjcZ+vVmhQf66+cLz7C4LqDLnNDDQCVjRAeCNUr7ZrGc+32CeP9G/ra7v0sTukgD3GhgIAPBcoy9urlsvOraS8/CsH/X+sky7SwKqhMu3rlw9mXzLli3nUg8AoJrc17uVDheWaNLCbbpv2iozY2dIV1Z24KVBJz09XU2aNNHgwYPVoEGDqq0KAFDlfHx89M8/tzGdWBPmp5uVncKSMv21Z1O7SwOqP+hMnTpV48eP13PPPae+fftq5MiRuuyyy+Try90vAPDksPPoFfEK9PfVa99t0f/7eK1Z2Sm/rQV4OpdTysCBAzVnzhzTWp6UlKS///3vio2N1QMPPKBNmzZVbZUAgCoNOw/0aa07/3hsRMh/PluvF7/axJwdOMIZL8c0atRIDz30kAk3U6ZMMedetW7dWgcPHqyaCgEA1RJ27r60pe7tfWyuzvNfbdR/v9hA2IHHO+M5OpaCggJNmzbN3Mqygo612hMSElL51QEAqr0bK8jfV49/sk4p3/ykwuIyPXT5+SYIAY4POlaoGTdunN5//33ThWXt05k+fTpnXAGAg9x4QTOzZ8caKvjmD1vNBmVr07KvL2EHDg46bdq0UVZWlum6+u6775SQkFC1lQEAbDO0W5wC/Hz14MzVenvRNhWXlumJ/u1MhxbgyMnIVndVaGio/P39f3MJ88CBA/IkTEYGgN8+G+veaStlHXY+ILGRnr6mvfz96LaFA8+6mjBhQmXVBgDwoLOxrNtYd01doRnLd6iwtEwvDOpgVnsAT+By0Bk2bFjVVgIAcEv9EqJNsLnj3TR9smqXikvK9L/BiQry97O7NOB3uRTJndheaJ1cHh8fr+TkZLtLAQC316dtlF67Icms7nyxdo9ueTtVBcWldpcFVE7QsTYiv/feeyoqKvrN91mzdW699VY99dRTcnejR4/W2rVrtXTpUrtLAQCPcEnrSI0b1knBAb76ZsNe3fjWMh0pIuzAAZuR586dq/vvv98c2HnppZeqU6dOio6OVnBwsBkUaAWGH374QWvWrNHtt9+uBx980GwQ8gRsRgaAM7Noy36NnLhU+UWl6ty0jsYPT1bNoLMaywZU+c9vl7uuLFaYsc68mjdvnrZt26YjR46oXr16SkxMVO/evXX99dd73Ewdgg4AnLnUbQc0fPxS5RWWqGPjWpo4srPCgwPsLgteJLcqgo4TEXQA4OyszMzWDeMWK7egRO1jIjRpZGfVCgm0uyx4iVwXf37THwgAOCsJsbX07k1dVTskQKu25+i6NxZr/6FCu8sCTkDQAQCctTbREXrvpm6qVzNI63bl6trXFykrr8DusoAKBB0AwDlpFRWmqTd3VWR4kDZlHdK1ry3SPlZ24CYIOgCAc3Ze/Zp6/+ZualSrhrbsO6wRE5bqUGGJ3WUBBB0AQOVoUjdUb/+1s+qEBmr1jhzd/PYyFZYwZwceGHSKi4uVmZmpDRs2eNwhngCAqtOsfk1NGJ6skEA/zd+8X/e8v1Jl1omggLsHnby8PL3yyiu68MILTRtXXFyczj//fNWvX19NmjTRqFGjmDIMADDdWNZxEQF+Pvp41S79c/YaRx4lBAcFneeee84EG+sE8169emnWrFlasWKFNm7cqIULF+qxxx5TSUmJ/vSnP6lPnz7mKAgAgPe6oEV9/Xdggnn+1sJtSvlms90lwUu5NDDwuuuu08MPP2zOvPotBQUFmjhxogIDAzVy5Eh5AgYGAkDVmTB/q/45e615/tSAdrq2c2O7S4JDMBnZRQQdAKhaT3+2Xi9/+5N8faRXhiSpd5sou0uCA1TZZORvvvnmtB9LSUk5008HAHC4e3u30qBOsbL2JN/x7nIt3rLf7pLgRc446AwYMECpqaknXX/xxRc1ZsyYyqoLAOAQPj4+eqJ/W/U6P1JFJWW6cdIyM0UZcMug88wzz6hv375av359xbVnn31Wjz76qD755JPKrg8A4AD+fr56aXCikuNqK6+gRMPGL1HmgXy7y4IXOOOgc+ONN+of//iH6b5KT0/Xf/7zH/3rX//Sp59+qgsuuECewrrNFh8fr+TkZLtLAQCvEBzgpzeHJqtVZJiy8gpN2OEQUFS1s96MfP/992vcuHEqLS3VnDlz1LVrV3kiNiMDQPXanVOgq19ZoB3ZR9Q+JkLvjuqq0CB/u8uCQ39+u/Qna+zYsSdda9SokUJCQvSHP/xBS5YsMQ/LnXfeeS51AwAcLioiWJP+2lnXvLJAq7bn6JbJqRo3LFmB/pxKBJtWdJo2beraJ/Px0ZYtW+RJWNEBAHusyMzW4DcWKb+oVH9OiNYLgzrI1+pBB6p7RWfr1q2uvA0AAJd1iK1l5ur8deJSfbRyp+rWDNSjV8Sb/2kGKgvrhAAA21zY8pejIibMT9cr3/1kd0nwxqDz1FNPKT/ftTbAxYsX02YOAHDZVYmN9MgV8eb5059t0PtLM+0uCd4WdNauXWtOKL/ttttMh9XevXsrPmYd5rlq1Sq9/PLL6t69uwYNGqSwsLCqrBkA4DB/7dlUt1x4nnn+wIxV+nLtHrtLgjcFnUmTJumrr75ScXGxBg8erKioKHNwpxVogoKClJiYqPHjx2vo0KFmkKDViQUAwJm4v08rDUyKMUdF3D4lTUvTD9hdErxxjk5ZWZlWrlypjIwMHTlyRPXq1VOHDh3MPz0RXVcA4D5KSstMu/lX67IUHuyvD27prlZR3CXAyTi93EUEHQBwL0eKSjVk3GKlbjuoyPAgTb+1u2Jqh9hdFpx+erk1Adk67qFHjx7m2IQHHnjArOgAAFCZagT6adywTmoZWVN7cgs1dPwSHThcZHdZ8FAuB50nn3xSDz74oGrWrGmmIlunlY8ePbpqqwMAeKVaIYF6a2RnRUcEa8vewxo5canyi0rsLgtODjrWhmSrs+rzzz/XrFmzNHv2bL3zzjtmzw4AAJWtYUQNc1RErZAAM0X51slpKi7lZw6qKOhYm48vu+yyitfW6eXW9MqdO3ee4ZcEAMA1zRuEafzwZAUH+Oq7jXt137RVKrPasoDKDjrWvJzg4OATrgUEBJiWcwAAqkrHxrX1yvVJ8vP10czlO/TvOevsLgkexKWzrixWc9bw4cPN3JxyBQUFuuWWWxQaGlpxbcaMGZVfJQDAq13cuoGevrq97vlgpd6Yt1X1w4J00x+ODRgEKiXoDBs27KRrQ4YMcfWXAwBwTq5OitH+w4V68tP15lE3NMhcAyol6EyYMMHVtwIAUCWsVZy9eYVmVee+6atUJzTQrPYAp8Pp5QAAjzKm7/nqn9hIpWVHdds7aUrLOGh3SXBjBB0AgEfx9fXR09e010Wt6utIcamZsbM5K8/usuCmCDoAAI8T4Oerl6/vqITYWsrOL9bQcUu0K4dp/TiZ1wadlJQUxcfHm+MsAACeJyTQXxOGJ6tZ/VDtzCkwYSc7n6MicCIO9eRQTwDwaNsP5uuaVxZqd26BkprU1uS/djHnZcHZKv1QTwAA3JF1srl1LlZ4sL858fz2KWkq4agI/IygAwDweK2iwjRueLKC/H01d32WHpy52gy6BQg6AABHSI6ro5cGd5Svj/T+su16/qtNdpcEN0DQAQA4xqXxkXr8qnbm+di5m/Tukgy7S4LNCDoAAEcZ3KWx7rykuXn+8Kwf9fX6PXaXBBsRdAAAjvP3S1tqYFKMmZ48+p3lWpGZbXdJsAlBBwDgOD4+PnpyQDtd2PKX6cnp+w7bXRZsQNABADh6enLbRuE6cLhIwyYs0b5DhXaXhWpG0AEAOFZokL/GD09WTO0a2rY/X3+duFT5RSV2l4VqRNABADhag7BgM1CwdkiAVm7P0e1TljNQ0IsQdAAAjnde/Zp6c9ixgYJfr8/SIx/+yEBBL0HQAQB4BescrLHXJZqBgu8uydT/vt5sd0moBgQdAIDX6N0mSv+8sq15/tyXG/X+sky7S0IVI+gAALzKDV2b6LaLzjPPx8xYrW83ZNldEqoQQQcA4HXu7d1KAxIbmYGCt72TptXbc+wuCVWEoAMA8MqBgk9d3V49m9dTflGpRkxcooz9+XaXhSpA0AEAeKVAf1+9MqSj4huGa9+hYwMFrcGCcBaCDgDAa4UFB2jCiGQ1qlVDW/cd1qhJy1RQXGp3WahEBB0AgFeLDLcGCiYrPNhfqdsO6r5pq5ix4yAEHQCA12veIEyvDkmSv6+PPlq5U89/tcnuklBJCDoAAEjq3ryenuzfzjwfO3eTZqRtt7skVAKCDgAAP/tLcqxu/XnGzv3TV2nxlv12l4RzRNABAOA49/6plS5v11DFpUd18+RUs0kZnougAwDAcXx9ffTsXxLUIbaWsvOLNWLCEh2k7dxjEXQAAPiV4AA/vTG0k2Jq11D6/nzd/HaqCktoO/dEBB0AAE6hfliQJgxPVliQv5akH9AD01fTdu6BCDoAAJxGi8gwvTyko/x8fTRz+Q6NnbvZ7pJwhgg6AAD8hgta1NfjV7U1z5//aqM+XLHD7pJwBgg6AAD8jus6N9bNf2hmnt/7wSotTT9gd0lwEUEHAAAX3N+ntfq0iVJRaZlumrRM6bSdewSvDTopKSmKj49XcnKy3aUAADyk7fz5QR2UEBOhg/nFGjlxqbLzaTt3dz5HvXwLeW5uriIiIpSTk6Pw8HC7ywEAuLmsvAL1T1mgHdlH1LVZHU0a2UWB/l67buD2P7/5nQEA4Aw0CAvW+OHJqhnkr0VbDmjMDNrO3RlBBwCAM9QqKkwp1x9rO5+etl0p39B27q4IOgAAnIULW9bXP//cxjz/7xe0nbsrgg4AAGdpSNcmGnVBU/P83mmrlLqNtnN3Q9ABAOAcPND3fP0pPlJFJWUaNSlV2/bTdu5OCDoAAJwDa5/OC9d2ULtGETpwuEgjJi5VTn6x3WXhZwQdAADOUUigv8YN66ToiGBt2XtYN09eZlZ4YD+CDgAAlaBBeLDGj/il7fzBmbSduwOCDgAAlaR1VLheGpxobmdNS6Xt3B0QdAAAqEQXtWqg/zuu7Xz2yp12l+TVCDoAAFSyG7o20Y09j7Wd3/PBStrObUTQAQCgCoy57HxdStu57Qg6AABUAWufzou0nduOoAMAQDW1nd/6TqqKS2k7r04EHQAAqrjtfNzwZIUG+mnBT/v1z9lr7C7JqxB0AACoYuc3DNeL1ybKx0eavChDkxam212S1yDoAABQDXrFR+r+Pq3N83/OXqt5m/baXZJXIOgAAFBNbv5DMw3o2EilZUc1+p00bdl7yO6SHI+gAwBANfHx8dG/B7RTx8a1lFtQohvfWkYnVhUj6AAAUI2C/P302g2d1KhWDW3Zd1ijp6TRiVWFCDoAAFSz+mFBemNoJ4UE+umHzfv0+Mdr7S7JsQg6AADYID46XM8P6mCev7VwmyYv2mZ3SY5E0AEAwCa920Tp3t6tzPPHPlqjBZv32V2S4xB0AACw0W0Xnaf+icc6sW59J01b93EmVmUi6AAA4AadWImNaynnSLH++tZS809UDoIOAAA2Cw6wOrGS1PDnM7Fun5KmEjqxKgVBBwAAN9AgLNh0YtUI8NO8Tfv0+Cfr7C7JEQg6AAC4ibaNIvT8oATzfOKCdE1ZnGF3SR6PoAMAgBvp07ah7rm0pXn+6Ic/auFP++0uyaMRdAAAcDO3X9Jc/RKiVWI6sVKVsT/f7pI8FkEHAAA37MR65pr2SoiJUHb+sU6svAI6sc4GQQcAADftxHp9aCdFhQdrU9Yh3fnucjNrB2eGoAMAgJuKDD/WiRUc4KtvNuzVU3PoxDpTBB0AANxYu5gI/XfgsU6sN+Zt1fvLMu0uyaMQdAAAcHNXtI/W3/7Ywjx/aOZqLU0/YHdJHoOgAwCAB7CCzuXtGqq49KhufjtVmQfoxHIFQQcAAA/g6+tjbmG1bRSuA4eLdONby3SosMTustweQQcAAA9RI9DPbE5uEBakDXvydNd7dGL9HoIOAAAepGFEDdN2Hujvq6/WZenZLzbYXZJbI+gAAOBhOsTWMgMFLS9/+5M+WrnT7pLcFkEHAAAPdGWHRrr5wmbm+X3TVurHHTl2l+SWCDoAAHio+3q31sWt6quguEyjJi3T3rxCu0tyOwQdAAA8lJ+vj168LlHN6odqV06BbpmcqsKSUrvLcisEHQAAPFh4cIDeHNpJYcH+St12UI/OWqOjR+nEKkfQAQDAwzWrX1P/uy5Rvj7S1GWZemtBut0luQ2CDgAADnBRqwZ6oG9r8/z/fbJO8zfvs7skt0DQAQDAIUZd0EwDEhuZIYKjp6Rp2/7D8nYEHQAAHMLHx0dPDminhNhays4vNp1Y3n5MhMcHnezsbHXq1EkdOnRQ27Zt9cYbb9hdEgAAtgkO8NPrNySZYyI27jmkv09doTIvPibC56iHb80uLS1VYWGhQkJCdPjwYRN2li1bprp167r063NzcxUREaGcnByFh4dXeb0AAFSH5RkHNej1RSoqKdMdlzTXPX9qJSdx9ee3x6/o+Pn5mZBjsQKPlds8PLsBAHDOEhvX1r/7tzPP//f1Zs1ZvUveyPag8/3336tfv36Kjo429xZnzZp10ntSUlIUFxen4OBgdenSRUuWLDnp9lVCQoJiYmJ07733ql69etX4HQAA4J6uTorRX3s2Nc/v+WCl1u3KlbexPehYt5uskGKFmVOZOnWq7r77bj322GNKS0sz7+3du7eysrIq3lOrVi2tXLlSW7du1ZQpU7Rnz57Tfj1r1cda7jr+AQCAU43p21o9m9dTflGpbnp7mQ4eLpI3sT3o9O3bV48//rj69+9/yo8/99xzGjVqlEaMGKH4+Hi9+uqr5lbV+PHjT3pvZGSkCULz5s077df797//be7plT9iY2Mr9fsBAMCd+Pv5mmGCjeuEKPPAEdN2XlJaJm9he9D5LUVFRUpNTVWvXr0qrvn6+prXCxcuNK+t1Zu8vDzz3NqQZN0Ka9Xq9BuuxowZY95X/sjMzKyG7wQAAPvUDg3UG0M7KSTQTwt+2q8nPl0nb+HWQWffvn2mq8paqTme9Xr37t3m+bZt23TBBReYlRzrn3fccYfatTu2+epUgoKCzO7s4x8AADhdq6gwPfeXBPN8wvx0fbDMO/5H318ernPnzlqxYoXdZQAA4Pb6tG2oO//YQmPnbtJDM39U8wY1TXeWk7n1io7VPWW1j/96c7H1Oioqyra6AADwVHf9sYUujY9UUWmZbn47VXtyC+Rkbh10AgMDlZSUpLlz51ZcKysrM6+7detma20AAHgiX18fPT+og1pG1lRWXqFumZyqwpJSOZXtQefQoUPm1lP57SerRdx6npGRYV5breXWsQ5vvfWW1q1bp1tvvdW0pFtdWAAA4MzVDPLX6zd0Uniwv5ZnZOvhmT86dtiu7Xt0rOMaLr744orXVrCxDBs2TBMnTtSgQYO0d+9ePfroo2YDsnWm1WeffXbSBmUAAOC6uHqhemlwRw2fsEQfpG5Xu5gIDe0WJ6fx+LOuzhVnXQEAvNkb328x7eb+vj6afGMXdW3m2lmRdvOas67OljWJ2RpAmJycbHcpAADY5sYLmurKDtEqKTuq0e+kaUf2ETkJKzqs6AAAvNyRolJd/coCrd2Vq7aNwjXtlu4KDvCTO2NFBwAAuKRGoJ9eH5qkOqGB+nFHrsbMWO2YzckEHQAAoJjaIXppcKL8fH00c/kOjfthq5yAoAMAAIzu59XTw5efb57/e856zd+8T56OoAMAACoM7x6nAR0bqbTsqG6fkqbMA/nyZAQdAABQwcfHR0/2b6f2MRE6mF+sm95OVX5RiTwVQQcAAJzA6rh6dUiS6tUM1Lpdubpv2iqP3ZxM0AEAACeJrlVDL1+fZAYJfrxql17/fos8kdcGHQYGAgDw2zo3raPH+sWb5//5bL3mbdorT8PAQAYGAgBwWlZMsG5dWedh1QoJ0Ozbeyq2TojsxsBAAABQKZuT/99Vbc3m5OyfNydbk5Q9BUEHAAC4tDm5buixzcljZnjO5mSCDgAAcGlzcsr1Hc3k5Fkrdmr8/HR5AoIOAABwSddmdSsmJz/56Tot+Mn9JycTdAAAwJlNTk4sn5y8XDuyj8idEXQAAMCZTU4e0E5tosN14HCRbnk7VQXF7rs5maADAADOeHPyazckqXZIgFbvyNGDM1e77eZkgg4AADhjMbVDlDK4o3x9pBlpOzRp4Ta5I4IOAAA4K92b19OYvsc2J/+/j9dqydYDcjdeG3Q4AgIAgHN34wVN1S8hWiVlR3XbO2nanVMgd8IREBwBAQDAOckvKtGAlxdo/e48JTaupfdu6qogfz9VJY6AAAAA1SIk0N9sTg4P9tfyjGz9c/ZauQuCDgAAOGdN6obqxesS5eMjTVmcoalLM+QOCDoAAKBSXNyqge65tKV5/sisNVqRmS27EXQAAEClue2i5vpTfKSKSst06+RU7TtUKDsRdAAAQKXx9fXRs39J0Hn1Q7Urp0Cj30lTcWmZffXY9pUBAIAjhQUH6LUbOqlmkL8Wbz2g7zbsta0Wf9u+MgAAcKzmDWrq+UEdVFRSpl7xkbbVQdABAABV4lIbA045bl0BAADHIugAAADH8tqgw1lXAAA4H2ddcdYVAAAeh7OuAACA1yPoAAAAxyLoAAAAxyLoAAAAxyLoAAAAxyLoAAAAxyLoAAAAxyLoAAAAxyLoAAAAx/L608vLB0NbExYBAIBnKP+5/XsHPHh90MnLyzP/jI2NtbsUAABwFj/HraMgTsfrz7oqKyvTzp07FRYWJh8fH1tqsA4WXbp0qZzK3b8/u+qrrq9bFV+nsj5nZXyes/kc1v8JWv9zk5mZyRl3NnH3/y44/ftLtrG+yvraVnyxQk50dLR8fU+/E8frV3SsfzkxMTG21uDn5+fo/9i6+/dnV33V9XWr4utU1uesjM9zLp/D+nXu/GfTydz9vwtO//78bKyvMr/2b63klGMzshsYPXq0nMzdvz+76quur1sVX6eyPmdlfB53//MF7/x9c/fvb7SN9VX31/b6W1cAvI9168r6P8GcnBy3/r9uAOeOFR0AXicoKEiPPfaY+ScAZ2NFBwAAOBYrOgAAwLEIOgAAwLEIOgAAwLEIOgAAwLEIOgAAwLEIOgDwO/r376/atWvrmmuusbsUAGeIoAMAv+Nvf/ubJk2aZHcZAM4CQQcAfsdFF11kDv4F4HkIOgA82vfff69+/fqZE4x9fHw0a9ask96TkpKiuLg4BQcHq0uXLlqyZIkttQKofgQdAB7t8OHDSkhIMGHmVKZOnaq7777bHPmQlpZm3tu7d29lZWVVvKdDhw5q27btSY+dO3dW43cCoCpwBAQAx7BWdGbOnKmrrrqq4pq1gpOcnKyXXnrJvC4rK1NsbKzuuOMOPfDAAy5/7m+//dZ8jmnTplVJ7QCqBis6AByrqKhIqamp6tWrV8U1X19f83rhwoW21gagehB0ADjWvn37VFpaqsjIyBOuW693797t8uexgtHAgQP16aefKiYmhpAEeBB/uwsAAHf31Vdf2V0CgLPEig4Ax6pXr578/Py0Z8+eE65br6OiomyrC0D1IegAcKzAwEAlJSVp7ty5FdeszcjW627dutlaG4Dqwa0rAB7t0KFD2rx5c8XrrVu3asWKFapTp44aN25sWsuHDRumTp06qXPnznrhhRdMS/qIESNsrRtA9aC9HIBHs9q+L7744pOuW+Fm4sSJ5rnVFv7MM8+YDcjWzJyxY8eatnMAzkfQAQAAjsUeHQAA4FgEHQAA4FgEHQAA4FgEHQAA4FgEHQAA4FgEHQAA4FgEHQAA4FgEHQAA4FgEHQBV6qKLLtJdd911xr/ukUce0U033SRP8+qrr6pfv352lwHgZwQdAG7HOqrhxRdf1EMPPVRxbfjw4fLx8dEtt9xy0vtHjx5tPma9pyoVFBSYr9GuXTv5+/vrqquuOuk9I0eOVFpamubNm1eltQBwDUEHgNt588031b17dzVp0uSE67GxsXrvvfd05MiRE8LHlClTzAGeVa20tFQ1atTQnXfeqV69ep32xPTBgweb87QA2I+gA6BaffLJJ4qIiNA777xz2vdYYeZUt386duxows6MGTMqrlnPrZCTmJh40i2z22+/3Tysr1evXj1zO+z44/0KCwt1//33m88ZFBSk5s2ba9y4caetKzQ0VK+88opGjRqlqKio077Pqv2jjz46IZABsAdBB0C1sVZerrvuOhNyrr/++lO+58CBA1q7dq06dep0yo9bt4YmTJhQ8Xr8+PEaMWLEKd/71ltvmVtMS5YsMbfCnnvuObNaVG7o0KF69913zerLunXr9Nprr6lmzZrn/H1atZeUlGjx4sXn/LkAnBv/c/z1AOCSlJQUs+dm9uzZuvDCC0/7voyMDLPqEh0dfcqPDxkyRGPGjNG2bdvM6/nz55sVoG+//fak91orNc8//7zZv9OqVSutXr3avLZWZDZu3Kj3339fX375ZcVtqGbNmlXK9xoSEmJWkcprBGAfgg6AKjdt2jRlZWWZUJKcnPyb7y2/3RMcHHzKj9evX1+XX365Jk6caAKR9dy6LXUqXbt2NSGnXLdu3fTss8+avTYrVqyQn5/faUNXmzZtKoLKBRdcoDlz5uhMWHt58vPzz+jXAKh8BB0AVc7aP2N1Ilm3mazbOseHj18rDy0HDx40oeZ0t6+svTflK0Vnwwoiv+XTTz9VcXGxS+893S2409UPoPqwRwdAlTvvvPP0zTff6MMPP9Qdd9zxu+8NDw83+3ROp0+fPioqKjJBpHfv3qd936/3yCxatEgtWrQwKzlWi3hZWZm+++67U/5aq+PL2pxsPRo1aqQz8dNPP5lusF9vkAZQ/Qg6AKpFy5YtTdiZPn36bw4Q9PX1NXtmfvjhh9O+xwoq1uZhKwxZz39rv8/dd9+tDRs2mE3H//vf//S3v/3NfCwuLk7Dhg0zq0OzZs3S1q1bzT4fa9/Ob7G+pnXby1qxycnJMc+tx/GsGTrWfh8rtAGwF7euAFQba0Pw119/bVq/rYBi7Zc5lRtvvNFsGH766adN8DkVa9Xn91hdVdaen86dO5uvZ4Wc46ctW63iDz74oG677Tbt37/ftKlbr3/LZZdddsIm4/JVm+Pb1q1QZdUPwH4+R4//2wkAbsD6z1KXLl3097//3bSjnw0rTHXo0EEvvPCCqtOaNWt0ySWXmK4uq/MKgL24dQXA7ViblV9//XUzi8bT7Nq1S5MmTSLkAG6CW1cA3JK1GmM9PM3pjoYAYA9uXQEAAMfi1hUAAHAsgg4AAHAsgg4AAHAsgg4AAHAsgg4AAHAsgg4AAHAsgg4AAHAsgg4AAHAsgg4AAJBT/X8XlA99/CGIegAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(keff, power_model_1d_1)\n", "plt.xscale(\"log\")\n", "plt.yscale(\"log\")\n", "plt.xlabel(\"k (Mpc-1)\")\n", "plt.ylabel(\"P(k) (Mpc3)\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "0be505f1", "metadata": {}, "source": [ "Now let's see what will happen if you naively use the survey volume and do a cubic box with some uniform k-sampling:" ] }, { "cell_type": "code", "execution_count": 13, "id": "48531876", "metadata": {}, "outputs": [], "source": [ "ps_temp = PowerSpectrum(\n", " wproj=ps.wproj,\n", " num_pix_x=ps.num_pix_x,\n", " num_pix_y=ps.num_pix_y,\n", " nu=ps.nu,\n", " ra_range=ps.ra_range,\n", " dec_range=ps.dec_range,\n", ")\n", "# cubic box with the same volume\n", "ps_temp._box_len = np.ones(3) * np.prod(ps.box_len)**(1/3)\n", "# same random grid resolution\n", "ps_temp._box_ndim = np.array([400, 400, 400])\n", "ps_temp.propagate_field_k_to_model()\n", "\n", "ps_temp.k1dbins = np.linspace(0, 0.5, 21)\n", "power_model_1d_2, keff2, nmodes2 = ps_temp.get_1d_power(ps_temp.auto_power_matter_model)" ] }, { "cell_type": "code", "execution_count": 14, "id": "9d0a534f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAG1CAYAAADwRl5QAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWeNJREFUeJzt3QV4k+caBuCnqTstpUYFlwIrVhw2NhiuY7jr0G1MmQtTNgYMhkNxt+EDhmtpKa5DCpQKlLo3Pdf35RQZsgBt//zJc19XDn/SNHl7WNuHz16z3NzcXBAREREZIY3SBRAREREVFAYdIiIiMloMOkRERGS0GHSIiIjIaDHoEBERkdFi0CEiIiKjxaBDRERERotBh4iIiIyWBUycVqtFZGQkHB0dYWZmpnQ5REREpAdx3nFSUhK8vb2h0Tx53Mbkg44IOb6+vkqXQURERM/h+vXr8PHxeeLHTT7oiJGcvP+jnJyclC6HiIiI9JCYmCgHKvJ+jz+JyQedvOkqEXIYdIiIiNTlv5adcDEyERERGS0GHSIiIjJaJj91RURExicnJwdZWVlKl0EvwNLSEubm5nhRDDpERGRUW46joqIQHx+vdCmUD4oUKQJPT88XOv6FQYeIiIxGXshxd3eHnZ0dz0dTcWBNTU1FTEyMvO/l5fXcr2WyQWfKlCnyJoY3iYhI/cTP87yQU7RoUaXLoRdka2sr/xRhR/ydPu80lskuRh4+fDjOnDmDkJAQpUshIqJ8kLcmR4zkkHGw+//f5YustzLZoENERMaJ01XGwywf/i4ZdIiIiMhoMegQERGR0WLQISIiIpQoUQITJkyAsWHQISIiUtgrr7yCd955R+kyjBKDTgHZcyEWJ06EIfdoMJCZonQ5RERkBGfLZGdnK12G6jDoFNB/jN9sOIPQ5T/CbMPbyPi5PLI2jQHiLitdGhGR6R08l5mtyE28tz769u2L3bt3Y+LEiXKXkbhdvXoVu3btktebN29GjRo1YG1tjX379snnt2/f/qHXEKNBYlQoj1arxQ8//ICSJUvK82gCAwOxcuXK/6wlKSkJ3bp1g729PYoXLy7Pm3tQREQE2rVrBwcHBzg5OaFz586Ijo6WHzt37pzcDr548eJ7z1++fLl8f3Gci1JM9sDAgpSamYOa/i64Gu+Da1p3+GfHAEf+QO6RqUj1awz7hsOB0q8CGuZMIqKClJaVg4Avtiry3me+aQY7q//+NSsCzoULF1C5cmV888038rFixYrJsCN8/PHH+OWXX1CqVCm4uLjo9d4i5CxcuBDTpk1D2bJlsWfPHvTs2VO+7ssvv/zEzxs3bhw++eQTfP3119i6dSvefvttlCtXDk2bNpXhKS/kiGAmRpfEmXRdunSRoaxChQqyzmHDhqFBgwbQaDR466238NNPPyEgIABKYdApAPbWFvjxjZeQ0OInrDw6FBf3r0Xz1PV4xfw47CP+Bhb9jXjXQDgM2wkLixdvWEZEROrl7OwMKysrORoi+jr9mwg/ImjoKyMjA99//z22b9+OunXrysdESBKjQdOnT39q0Klfv74MVoIIOPv378dvv/0m33/Hjh04efIkrly5Al9fX/mc+fPno1KlSvLw3aCgIBlyNm3aJEOV+JrEYyNHjoSSGHQKkLOdJQY0Kgttg/ex71JffLp3P0pfWYpO5ruwNMYPwT/vQvfafugaVBzu6RGAewWlSyYiMiq2luZyZEWp984PNWvWfKbnX7p0SfaJ+nc4yszMRLVq1Z76uXnB6MH7eTuxzp49KwNOXsgRxEiNaLwpPiZCjTBnzhwZksSIzunTpxU/wJFBpxBoNGZoVK4YGpVrjxt3X8ecg+ew9ug1RCWmY/y2Cwj9exXmWf6ARI/acGw0DGYVWgHmlkqXTUSkeuKXrD7TR4ZMrJd5kAgQ/17/82CLhOTkZPnnxo0b5TqbB4l1PgXt+PHjSElJkXXeunXrhRpy5gd1/+2rkI+LHd5pWR1DXw/E5pNRWHDoGircuIrsXA2cog8DKw4jxdodlrUHwKrWAMChmNIlExFRARPTPPo2mRbrbE6dOvXQY+Hh4bC0tLw3yiICjVg4/PJTpqke59ChQ4/cr1ixorwWf16/fl3e8kZ1xCJj0Ug1bw1OXFycXCz96aefypDTo0cPhIWF3WvQqQQGHYVYW5ijfbXi8nY6shJ+2fMGipxZhE5mO+CWEQPs+QHZe39Batk2cHpjEmDtqHTJRERUgIf1HT58WC5AFot9XV1dn/jcV199VS4aFutjxNSSWHQsgk/etJSjoyPef/99vPvuu3IBsVgYnJCQINfbiJ1Sffr0eeJri+f8/PPPclfXtm3bsGLFCjkyJDRp0gRVqlSR4UVMZ4nFyGJNjghTedNrYvGxCEGfffaZXCskahK1/Hv3VmHith8DUMnbGR93bYpuY2Ziw2vb8b31uzimLQOL3CxEnTuMbsGnsPnkLWTlaMVeSaXLJSKifCbCgLm5uRwZESM2YjTmSZo1a4bPP/8cH374oVwXI7aE9+7d+6HnfPvtt/I5P/zwgxyJad68uQwsYrv507z33ns4evSoDChjx47F+PHj5fvlTQOuW7dO7vxq1KiRDD5ikfOyZcvkx0XwEguRFyxYAAsLCznlJkLYzJkz5RZ5pZjl6rvR30glJibKFe8i7Yqkawi02lwc+OcOdu/cgjNXI7FfW0k+7u+Yi7XmY2BVpR3s6w8GivgpXSoRkcFIT0+XO4LEL3MbGxuly6EC/jvV9/e3yU5diWE0cdN3TrSwFy83KOuGBmV74mZ8GpYcjsDSkAgEpW6Di2UEEPI7tCFTEO/bBC6vDIdZqZdF1Fa6bCIiIoPDER0DHNF5nIzsHGw9cQNnd69A/bjVaGB++t7H4u1Lwbb+EFjX6AlYOyhaJxGRUjiiY3zS82FEh2t0VLR4uW11f3z07vtwHboFv5VbgMXa15GSa40iKZdhufVjTFi3Dxejk5QulYiIyGCY7NSVmgV4OyGge1skpLXAmsNncffgfDilXMWEsBxMCNuDOqVc8bXLFpSuFASLCi0ADU9fJiIi08Sgo2LOtpbo+cpLyH15nFy83OzgVWw7E42IyxdQxnoizE/nItHaC2a1BsKxbn/A7snbFYmIiIwRg44REFv+6pdxk7fI+DSs23cMi0Lboq12O4pk3AL2fovMfT8hvnQ7FHt1BMy8qypdMhERUaHgGh0j413EFkNb10PXT4Oxv+1e/OH0Dk5r/WGVmwn3SytgNuNl7F09FckZ2UqXSkREVOA4omOkrCw0aFWjNFDja5y5ORrTdm6Ez8WFqI/jGHrEDbnh29Gxug8GlkmEv19JwEnZXiREREQFgSM6JiCguDPe6tkdjcasw4bXtsGjWDGkZOZgwaGrSFo+FNnjKyFyZldkX9nHk5eJiFQoODhYdhF/mq+++gpVqz7/0gXRw0q0hlAbBh0T4mRjiV6NArB99MtYPLA2OlRwRCpsYIEceN/cDIt5rRA7LggJ+2YCmSlKl0tERPncZmLHjh0wNQw6Jrp4uV4ZN/zW92X4vrcL8wMXYo1ZE6TlWqFY6kU4b38fqT+Wwz9b/oCJnydJRGQ0HBwcULRoUZgaBh0T5+Vsi94d2qDVpyuwu/VuzHcchAhtMdhpk/H9ntt4/bc9mH/wKpKSEkUTLqXLJSIySqLLuOgaXqZMGVhbW8PPzw/fffed/NiuXbvkP1Dj4+PvPT88PFw+JrqdP2jt2rUoW7asPEVYNOO8fv36U6eu5syZg0qVKsn39PLywogRI/6z1q+//lo2HhWnEYtu5ZmZmfc+JjqWjxo1Cu7u7rIG0Tk9JCTk3inH4r0GDx587/n//POP7LYu6igoXIxM9xYvNw8KAIJ+wbnITxG8bSUO/1McyTHJ+GLdaSRvWYXO1oehrTkQ7g37AbZPnwsmIjIoT5uONzMHLG30fK4GsLT97+da2T9TeWPGjJFdvn/77TcZDm7duoVz584902ukpqbKcCS6iFtZWWHYsGHo2rUr9u/f/9jnT506FaNHj8aPP/6IFi1ayFYKT3puHjH1JQKMCF8iZPXr10+OEuWFMtFRfdWqVZg3bx78/f1leBOB69KlS3B1dcWiRYtQu3ZttGrVCq1bt0bPnj3RtGlT9O/fHwWFQYceUcHbBRX6DMIb6VlYHXYTCw5eRbOEA3DLvAUc+ArpB39AlH97FH99JCy9qyhdLhHRf/ve+8kfK/s60GPF/fvjygBZqY9/rn8DoN/G+/cnVAFS7zz6vK8S9C4tKSkJEydOxOTJk9GnTx/5WOnSpWXgeRZZWVnyNUSQEETYqFixIo4cOYJatWo98vyxY8fivffew9tvv33vsaCgoKe+hwhQYvTFzs5Ojs588803+OCDD/Dtt98iLS1NhiexMFoEJ0GEt23btmH27NnyeWJESbzvwIEDZQi7du0aNmzYgILEqSt6IkcbS/SpVwLbRr+M2O5bsajYuziv9YVNbgZKXF0GyxkNcH18Y9wNW6t0qUREqnX27Fk55fPaa6+90OtYWFg8FFQqVKggd2KJ1/+3mJgYREZGPvN7BgYGypCTp27dukhOTpZTZGIaSoSt+vXr3/u4paWlDFkP1iDCVbly5WQoE6GpoNcNcUSH/pOYB65TwR91KnyFqPiPsHz7WhQ9HYyXtUfgmxiGJWvmY88ZX/Sq44+6pYvK5xMRGZRPIp8+dfWgDy495bn/Gh945+QLFgbY2j4wFfYYGo3uPR/cHCICRUG+Z0ESIevChQswNzfHxYsX0bx58wJ9P47o0DPxLGKLzp26odFnm7GrxQ6scuiGudmvY/OpKHSfdRjDxs3CP9O7I+XyIZ7JQ0SGQ6yZedLtwfU5//ncfwWEJz3vGYjFwyJ4PGnrt1j4K4h1Ow8uRv637OxsHD169N798+fPywXMYvrq38QC4BIlSjzzdvPjx4/LKao8hw4dkru5fH195XSbmNp6cJ2PCGRiMXJAQMC9x8R6nCpVqsiptY8++uixI075iSM69FwszTVoUqc6UGcaKkclycMH14TdRJOkdSiduheYvxE3bCtAU2cwvOv1ePQHCRERSWJxr/iFLxbyiqAgpn5iY2Nx+vRpDBgwQO7EEkFC7JoSi37FaMivv/76yOuIaaKRI0di0qRJchpL7KCqU6fOY9fnCOL1xK4psUNKrKkRa4VESBGv8SRih5Wo6bPPPpOLkb/88kv5PmLUyd7eHkOHDpVrccTCY7FzTCxGFoukxecIU6ZMwcGDB3HixAn5NW3cuBE9evSQgUl87QUi18QlJCSIYQf5J72YxLTM3A2bN+RuHds+N/2Lorm5XzrJW8JXvrkXFr2XmxF7RekSiciIpaWl5Z45c0b+qTY5OTm5Y8eOzfX398+1tLTM9fPzy/3+++/vfXzfvn25VapUybWxsclt2LBh7ooVK+TvritXdD9X586dm+vs7Jy7atWq3FKlSuVaW1vnNmnSJPfatWv3XuPLL7/MDQwMfOh9p02bllu+fHn5nl5eXrkjR458Yo19+vTJbdeuXe4XX3yRW7Ro0VwHB4fcQYMG5aanp997jvj/XryGm5ubrKF+/fq5R44ckR87e/Zsrq2tbe7ixYvvPf/u3bu5vr6+uR9++OEz/53q+/vbTPwPTFhiYiKcnZ3ltjpxJgC9OPGf1NEzF3Hz7xmodXs1vM10OxJuwB3L6q5H9zr+8vweIqL8JM5puXLlCkqWLClHSUj9nvZ3qu/vb05dUb4Ti5GDKpVDUKVfEB3/Ff7csghe5xfg78xKmLrzH/yx+zJer+CK94odRenGvWFmw4BJREQFw2SDjpgnFLecnBylSzFqHkUc0LbrEGTlDELMqSjUOXwNhy7HwezcRpS5PAkph8fimm97+DUfBYfi9xerERER5QdOXXHqqtBdiE5C6Jb5qHP5d5Q0u7+L4KJDEGzqD4Vv7faA5l/bPYmI/gOnroxPej5MXXF7ORW6ch6O6NZnOIqNOYHtNabhgEUtaHPNUDY5BL5b+yN6bEVsOnIamdnsrUVERC/GZKeuSHkONlZo0qYbclt3RfiJ44jdNRVBcRtwPdsZw1Zfhdtfkega5IeeARbw9C2tdLlERKRCDDpkEIuXqwVWBQKnI/rOXZw4eAweJ7SITszAwp3HMOLACFyyKYfsmgNRvnEPmFlYK10yERkwE1+RYVRy8+HvklNXZFA8irqgf+tXse+jV/FHj+ro4R0FDbQok3EaFfa/i7vflcOxeR8gMSZC6VKJyMCIA/MEcUAdGYfU//9d5v3dPg8uRuZiZIN3+fIlXP3rD1S+tRruZnflY1m55jhT5GXYtfgWZStUVrpEIjIQok2CaHsgTvsVzSfZe0+dRDQRIUf0xRKNSb28vJ779zeDDoOOaiSnpiFs6wK4nJqLKjlnkJlrjvoZv8PPvyR61/VH88qesLbgbi0iUyZ+pUVFRcmwQ+onQo6np+djAyuDjp4YdNRH/Cd7KnQ/ToTsxJfXayBbq/tPeI7tRDh6loZvs1HwLPFoEzsiMh3ijLQX7fBNyhLTVaLD+ZMw6OiJQUfdYhLTsTTkOg4cOoClWaPkY2Kr+in7WtDUHoyABh2geco3ChERqRODjp4YdIyD+Jfb8Z0rYRE6C1Uzjt57/IaZF26U7YGKzYfC2dVN0RqJiCj/MOjoiUHH+Fy7cByR26egUvR6OJnpVux/qB0GTWA39Krrj0rezkqXSEREL4hBR08MOsYrJSkep7bMgvX5deic/B4yodueOMz9FF4p747Apt1hbcUzeYiI1IhBR08MOsZP/CcecvUuFhy6hi0nb2K7xbvw18QgCkVxwfdNlG42DMV9/JUuk4iIngGDjp4YdExLTNxdXF3zLcpcXwlXJMjH0nMtscy5P+wbjUCLKt6wt+aB4UREho5BR08MOqYpOyMNp3csgMOxWSiddV4+dkhbEZ9jGF6qHIhONXxQu6QrNBoeNkZEZIgYdPTEoGPicnNxd+8M2O/6ElbaNHTP/AQHtLqTln1cbNGxug/eqF4c/kXtla6UiIgewKCjJwYdkuKuIPfiNoR5dsLK0BvYcPwW0jPSkfX/vre1SrjKUZ6WL3nBgVNbRESKY9DRE4MOPU5G7BXkzGmO+ba98dOtQOTm6qawbC3N0aKyJ96o4YO6pYpyaouISCEMOnpi0KHH2vQhcGS6vEwv3RxLPEZjwck0XI5NufeU4kXE1FZxvFHdByXcOLVFRFSYGHT0xKBDj5WTDeyfAOz6EdBmAbauyG39G8IdX5ZTW38ej0RSeva9p9f0d7k3teVkozuvh4iICg6Djp4YdOipok4Ca94Cok/p7lfuBLQch3RLZ2w/Gy1Dz54Lsfh/X1FYW2hkF3UxylO/jBvMObVFRFQgGHT0xKBD/yk7E9j9E7BvPJCrBV75BHjlo3sfjk5Mx9pjN2XouRiTfO9xTycb3dRWDR+ULuagUPFERMaJQUdPDDqktxuhuumsN2YBFo+2jhDfSiduJGBV2A2sC49EQlrWvY9V8ysip7Zav+QNZ1tObRERvSgGHT0x6NALreNZ2Q8IGgCUeuWhD2Vk52DH2RisCr2BXRdikfP/uS0rCw1eD/CQoadh2WKc2iIiek4MOnpi0KHndmgqsOVj3XXQIKDp14DVo7uvYpLSse5YpJzaOh+ddO9xd0drdKheHJ2q+6Csh2NhVk5EpHoMOnpi0KHnlpEMbPsCODpbd9+lJNBhGuBX57FPF99qpyMTZeBZG34T8an3p7YCfZzlKE+bQG8UsbMqrK+AiEi1GHT0xKBDL+zSDuDPkUDiTfEtBdQbCTT+FLC0eeKniKmtnedisDL0Jnaej7k/tWWuQZMAdxl6GpUtBgtzTSF+IURE6sGgoycGHcoXafHA1k+A8EW6+y91ATrO0OtTY5MysC5ct2vrXNT9qS03B2t0qOaNTjV8Ud6TU1tERA9i0NETgw7lq3ObgC0fAT1WAcXKPfOnn45MkIFH7NqKS8m893iV4s6yuWjbqsXhas+pLSKiRAYd/TDoUIHsxjJ/oPFnyCzAry7gUUnvl8jM1mLXeTG1dQN/n4tB9v+ntizNzfBaBQ95Ns8r5YvBklNbRGSiEhl0nm7KlCnylpOTgwsXLjDoUMG4fgSY0wwwMwcafwLUG/VwCNLDneQM2XJChB6xmDlPUXsrtK+m67UV4M3/donItCQy6OiHIzpUoJKigfVvAxc26+77BAHtpwFuZZ7r5c7eSpRn84hdW7eT709tBXg5yQXM7ap6o6jDo4cZEhEZGwYdPTHoUIET32LHlwCbPwIyEgELW6DJl0CtIYDm+aaesnK0sseWGOURPbeycnTfxhYaMzSuoNu11bi8uzygkIjIGDHo6IlBhwpNwg1g3Qjg8k7d/XLNgW5LAbMXOx35bkom1p/QTW2JFhR5xKLltoFi15YPKnk7wewF34eIyJAw6OiJQYcKlfh2OzoH+OtzoPn3QI2++fryF6KT5NTW6mM35bb1PBU8Hf8/tVUcxRw5tUVE6segoycGHVJEwk3Ayfv+aM7NUMDRS/dYPsjO0WLvxdtYGXYD205HIzNHKx8XvbUalXWTi5ibBnjAzurZFkYTERkKBh09MeiQ4tITgD/qAZlJQIufdYcN5uM0U3yqmNq6JUd6wq/H33vc1tIcr1fykAuYRYNRblUnIjVh0NETgw4pLj4CWN4HiAzT3a/QGmj9G+Dgnu9vdSkmGX+G38S645G4dif13uMudpZo9ZKXnNqq4ecCDbuqE5GBY9DRE4MOGcwhg/snALt+BLRZgF1RoNV4oFL7Ank78W0vRnfECcwbTkQ+tFW9eBFbtK3qLUd6Knjye4KIDBODjp4YdMigRJ0E1rwFRJ/S3a/cCWg/FbAouLYPYj3PgX/uyNCz9XQUkjOyH1rELEKP2L3l42JXYDUQET0rBh09MeiQwcnOBPb8DOwdD1RsA7wZnK9rdp4mPSsHO87GyCaju87H3lvELASVcJG9tlpV8WK/LSJSHIOOnhh0yGDdCAVcSgD2RXX3U+MAjQVgUzj/nSakZmHzqVtypOfQlTtyZ3zeoYSNyhWTU1vcuUVESmHQ0RODDqmC+DZd0Qe4GQa0mwyUeqVQ3/5WQho2HL+Fdcdv4tTN+/22uHOLiJTCoKMnBh1ShZQ7wMzGQPw13f2gQUDTrwEr+0Iv5VJMEv4Mj+TOLSJSFIOOnhh0SDUykoFtn+tOVhZcSgIdpgF+dRQpR5+dW+2rFkd5T0dF6iMi48agoycGHVKdSzuAP0cCiTfFtzBQbwTQ+DPA0kaxkrhzi4gKG4OOnhh0SJXS4oEtY4DjiwEnH2DYAcDGGYaAO7eIqDAw6OiJQYdU7dwmwNoBKNlId198O2uzAXNLGALu3CKigsKgoycGHTIqofOAkJlAh+mARyUYkqiEdKw/HvnEnVtiPU+Dsm7cuUVEemHQ0RODDhmNnCzg9+q63lkaS6DxJ0C9UYC54Y2W/NfOLRF6qnPnFhE9BYOOnhh0yKgkRQPr3wYubNbd9wkC2k8D3MrAEHHnFhE9LwYdPTHokNER39LHlwCbPwIyEgELW6DJl0CtIYDGcKeFuHOLiJ4Fg46eGHTIaCXcANYNBy7vAsw0wJA9gGcVqMGTdm6Jll8DG5TER80rwIJreYhMWiKDjn4YdMioiW/vo7OBlNvAKx9DjfJ2bq0Nv4lDl+PkY7VLumJy9+oo5mitdHlEpBAGHT0x6JDJib0A7PgaaDkOcPKGmmw+eQvvrziOlMwceDhZ448eNVDD30XpsojIgH9/c+yXyNRseBc4twH4ow5wfJlu1EclWlTxwroRDVDG3QHRiRnoOuMgFhy8Khc1ExE9DoMOkalpPR7wrg6kJwBrBgPLegLJsVALEXLWDq8vT1fOysnF5+tO473lx5GWmaN0aURkgBh0iExNsfLAgG3Aq5/rztsRoztT6wFX9kItHKwtMLl7NXzasiLMNWZYfewmOvyxH9fupChdGhEZGAYdIlMkDhFs9D4weCfgHgCkxADz2+p2aKmEmZkZBjUqhYUDasPNwQrnopLQ+vd92HE2WunSiMiAMOgQmTKx3XzgDiCwu+5wQf/6UJu6pYti/cgGqOZXBEnp2Rgw7yjGb7uAHC3X7RARd11x1xWRIH4MZKUCVvb320ncvgh4BEAtMrO1GLvxDOYfvCbvv1yuGCZ2rYoiduySTmSMuOuKiPQnTuLLCznC9q+AGS8DR+eoZleWlYUG37SrjPGdA2FjqcHuC7FyKuvUzQSlSyMiBTHoENHDtDnA3atATqZuK/qaIUCmehb5dqzug9VD68PP1Q437qbhjakHsDL0htJlEZFCGHSI6GEac6DLQqDpN4CZOXBiGTDzNd1BgyoR4O2E9SMa4NUK7sjI1spDBj9dcxIZ2dyCTmRqGHSI6PFTWfXfBvqsBxw8gNizwMzGwKlVUAtnO0vM6l0T7zYpJ7+cRYcj0Hn6IUTGpyldGhEVIgYdInqyEvWBIXuBEg2BzGTgz1GqOlxQozHD203KYk7fIDjbWuL49Xi5bufApdtKl0ZEhYRBh4ieztED6LUWaPge0HYS4FAMatO4vLucygrwckJcSiZ6zj6Mabv/YesIIhPA7eXcXk70fK7uB7LSgLJNoBbpWTn4dM0prArTLU5uXskT4958CY42lkqXRkTPiNvLiajgJEUDK/oCizoBO7/X7dRSARtLc/zy5kv4rkNlWJqbYcvpKLSbsh8Xo5OULo2ICgiDDhE9OxtnoGJrcdIgsPsnYGFHIOW2alpH9Kjtj+VD6sLL2QaXY1Nk2Nl44pbSpRFRAWDQIaJnZ2kDtP4N6DADsLTT9cia1hCIOAy1qObnIltH1C1VFKmZORi+OAzfbTyD7Byt0qURUT5i0CGi5xfYBRj0N1C0LJAUCQS3BA5OUc1pym4O1lgwoBaGvFxK3p+59wp6zDqM2KQMpUsjonzCoENEL8a9oq4LeqWOgDYbuBkGNbEw12BMi4qY2qM67K3McfhKHFr/vheh1+KULo2I8gF3XXHXFVH+ED9Kji8BKrYFrB2gRpdikvHWwlD5p1is/HnrAPSq4y/X9RCRYeGuKyIqXCIMVO1+P+SI4LNqIBC+GGpRxt0Ba4fXR6sqXsjKycUX605j9PLjSMtUx64yInoUgw4RFYzTa4CTK4C1Q4E/RwJZ6VADB2sLTO5eDZ+1qghzjRnWHLuJDn/sx7U76mlsSkT3MegQUcEIaA80/lQM9QBh84HZTYG4y1ADMVU1sGEpLBpYG24OVjgXlSRbR+w4G610aURkqkEnNTUV/v7+eP/995UuhYgEjQZ4+UOg1xrArigQdQKY/gpwbiPUok6potgwsiGq+xVBUno2Bsw7ivF/nUeO1qSXNhKpitEEne+++w516tRRugwi+rfSjXWNQX1rAxkJwNLuwL4JUAtPZxssHVwXvev6y/uT/r6E/sEhiE/NVLo0IjKVoHPx4kWcO3cOLVq0ULoUInoc5+JA341A3RGAxhLwrw81sbLQ4Jt2lTG+cyBsLDXYfSFWTmWduBGvdGlEZOhBZ8+ePWjTpg28vb3lvPjatWsfec6UKVNQokQJ2NjYoHbt2jhy5MhDHxfTVT/88EMhVk1Ez8zcEmj2HTDyKOAbdP/xlDtQi47VfbB6aH34udrhxt00vDH1AIL3X2EXdCIDpnjQSUlJQWBgoAwzj7Ns2TKMHj0aX375JcLCwuRzmzVrhpiYGPnxdevWoVy5cvJGRCrgUuL+ddQpYOJLwJ5fAK06Wi8EeDth/YgGeD3AQ25B/2r9GQxdGIaEtCylSyMiQz8wUIzorFmzBu3bt7/3mBjBCQoKwuTJk+V9rVYLX19fjBw5Eh9//DHGjBmDhQsXwtzcHMnJycjKysJ7772HL7744rHvkZGRIW8PHjgkXo8HBhIp4O+xwJ5xuuuyzYAO0wA7V6iB+NEZfOAqvt90VgYeX1dbTO5WHYG+RZQujcgkJBrDgYGZmZkIDQ1FkyZN7j2m0Wjk/YMHD8r7Ysrq+vXruHr1Kn755RcMGjToiSEn7/ni/5i8mwg5RKQQsf287e+AuTVwcSsw42XVtJAQ/zDrV78kVr5VT4ac63Fp6DTtAObs41QWkSEx6KBz+/Zt5OTkwMPD46HHxf2oqKjnek0xAiTSX95NhCQiUvA05eq9gYHbdFNa8RHAnGZAyGzVNAYVIzhiC3rzSp5yZOebDWcwZEEoElI5lUVkCAw66Dyrvn37ylGdp7G2tpZDXA/eiEhhXoHA4N1A+VZATiawcTRw5tGNCYbK2dYSU3tWx9dtK8HKXIO/zkSj5aS9CL/OXVlESjPooOPm5ibX3kRHP3waqbjv6empWF1EVABsiwBdFwFNvwXKNNE1B1URMZXVp14JrBpaT+7Kuhmfhk5TD2DW3sucyiJSkEEHHSsrK9SoUQM7duy495hYjCzu161bV9HaiKiAprLqjwK6rwA05rrHRI+sC1uhFlV8nLFhVAO0rOKJbG0uxm48i0HzQ3nAIJGpBh2xUyo8PFzehCtXrsjriIgIeV9sLZ85cybmzZuHs2fPYujQoXJLer9+/RSunIgKtH1Enq1jgMWdgU0fANn3d0waMicbS0zpXh3fttNNZW0/G41Wk/YhLOKu0qURmRzFt5fv2rULjRs3fuTxPn36IDg4WF6LreXjxo2TC5CrVq2KSZMmyW3nhbk9jYgUIH48iS3oe/+/9q54DeDNYKCIH9Ti1M0EDF8chmt3UmGhMcNHzStgYMOScqqLiAr+97fiQUdpDDpEKiCmrlYPBtLjAVsXoOMsoOz9YycMXVJ6Fj5efRIbT9yS91+r4I5f3gyEi72V0qURqZZRnKNTkMRJzAEBAfIwQiIycOWaAUP2AN7VgLS7wKJOwP6JUAtHG0tM7lYNY9tXln2zdpyLQatJexF6jVNZRAWNIzoc0SFSD7FGZ8sY4OhswMoBGBECOHlDTU5HJmDE4mO4cjtFTmV90Kw8BjUsBY2GU1lEz4JTV3pi0CFSITGa41ML8Ffn7ksxlfXJmlNYfzxS3m9cvhh+7VwVrpzKItIbg46eGHSIjEBSFODgoduerhLiR++SI9fx1frTyMzWwsvZBr93q4aaJdTR64tIaVyjQ0SmIeok8EcdYOd3qmkbIYhdV91r+2HtsPoo5WaPWwnp6DLjEKbu+gdarXq+DiJDx6BDROoWcUi3QFl0Qd/xtarCjhDg7YQ/RzZAu6reyNHm4qct59B/XgjuJKvjzCAiQ8egQ0TqVmsQ0PxH3fW+34BtX6gu7DhYW2BCl6r4sWMVWFtosOt8rDxg8MiVOKVLI1I9Bh0iUr86Q4EW43TXByYBf32murAjprK61vLD2uH1UaqYPaIS09Ft5iFM2XmJU1lEL4BBh4iMQ+3BQKtfddcHJ+u2oass7AgVvZywfkQDdKhWXE5ljdt6Hn2DOZVF9LwsnuXJoqHm7t27sXfvXly7dg2pqakoVqwYqlWrhiZNmsDX1xdqOjBQ3HJycpQuhYjyS9BAwMwc2PAOEHVCd+6OpQ3Uxt7aAuM7B6JuqaL44s9T2HMhFi0n7cXCAbVR1sNR6fKIVEWv7eVpaWn49ddfMXXqVMTFxcl+U97e3rC1tZX3T506hcjISLz++uv44osvUKdOHagFt5cTGaFzm4CSjQBrB6jd+agkDFsUin9iU+DmYI2lg+ugjLv6vy4igzpHR4zU1K1bF3379kXTpk1haWn5yHPECM/ixYsxffp0fPrppxg0aBDUgEGHyASc3wyUbfZwV3QVuZuSKdfrnItKgrujLuyUKsawQ6YtMT+DztmzZ1GxYkW93jgrKwsREREoXbo01IBBh8jI7R4H7BwLVO8NtJ6o2rAj1uh0n3kY56OT4OFkjWWD66KEm73SZREZx4GB+oYcQYz2qCXkEJEJcPEHzDRA2HzgzxGAVp3r8oo6WGPRoNoo6+6A6MQMOcITcSdV6bKIDN4z/dNGDP5cuXIF2dnZ8n5mZiaWLVuG+fPn4/bt2wVVIxHR83upM9Bxpi7shC8C1g5TbdgRa3QWD6qD0sV0JymLsHM9jmGHKF+Czvnz51GyZEmUKVNGjvCIwFOvXj0MGDAAQ4cOlY9dvHhR35cjIio8VToBb8zW7cg6sRRY8xaQo/sHm9oUc7TGkkF1ZNuIm/Fp6D7rkPyTiF4w6Hz00UcIDAxEeHg4WrdujVatWsHHxwd3796VO6/EYuVvvvlG35cjIipclTsCneYAGgvg5HJg7VuqPGdHcHeykSM7/kXtcD0uDd1mHMKtBIYdohfqXu7u7o6//vpLbi1PSUmBo6Mj9uzZgwYNGsiPHzhwAN26dZO7r9SEi5GJTMzZ9cCKfkDzH3TtI1QsMj4NXWYclGGnRFE7LBtSFx5O6js3iMggupcnJyfD1dVVXtvb28ubl5fXQ1vQo6Ojn6tYIqJCU7ENMDJU9SFH8C5iK6exfFxscfVOqhzZiUlMV7osIoOid9ARBwSKbeN5fv75ZznKkyc2NhYuLi5QC3EqckBAAIKCgpQuhYiU2ImVJzVO1wg0OxNq5ONiJ8NO8SK2uHw7Bd1nHUZsEttFED1z0BEtHs6dO3fvvliALKav8ohprerVq0Mthg8fjjNnziAkJETpUohIKWLmfnEXYP9EYEVf1YYdX1c7LB5UG17ONrgUk4wesw6xNxbRs67R+S9iF5aNjc1D01lqwDU6RCbu0nZgSXcgJwMo1xzoPB+wsIYaXbmdgq4zDspzdip4OsoFy672VkqXRaSONTr/RWw9V1vIISJCmSZA92WAhQ1wYQuwrCeQpc51LiXd7GW4EVvQRbuIHrMOIz5VnaNURPnlmYPOqFGjMGnSpEcenzx5Mt555538qouIqPCUbgx0Xw5Y2AIX/wKW9VBt2CldzAFLBtWWhwuevZWInrMPIyE1S+myiNQTdFatWoX69es/8rg4PHDlypX5VRcRUeEq9TLQYwVgaaebzto4GmpVxl1MW9VGUXsrnLqZiF5zDiMhjWGHTNMzB507d+7IObF/E/NjbANBRKpWsiHQYyXgVg5o9D7UrJyHo+yN5WJniRM3EtBnzhEkpTPskOl55qAjWkBs2bLlkcc3b96MUqVK5VddRETKKFEfGHoQcH3g55lKT1Cu4OmERQProIidJcKvx6Pv3BAkZ6iz9QXR87J41k8YPXo0RowYIc/NefXVV+VjO3bswK+//ooJEyY8dyFERAbD/IEfjRe2Aof+ALosAqwdoDYB3k5YOKA2us88hNBrd9Fv7hEE96sFe+tn/vFPZDrby6dOnYrvvvsOkZGR8n6JEiXw1VdfoXfv3lAbbi8noifKTAEmBgIpsYBfXd0aHuv754epyYkb8XIXVlJ6NmqVdEVwvyDYWTHskHrp+/v7hc7REaM6tra2cHBQ379y8jDoENFT3QgFFnQAMhIA39q6NTw26vxZIaaveomwk5GNuqWKYk7fINhamStdFpFhnqMTExMjT0o+duyYDDxqwxYQRKQXnxpA77WAjTNw/TCwsCOQngA1qupbBMH9a8HeyhwHL9/BoPlHkZ6Vo3RZRAXqmUd0kpKSMGzYMCxZsgRarVY+Zm5uji5dusjw8LgdWYaMIzpEpJfIcGB+OyA9HiheA+i5GrAtAjUKuRond2GlZuagUblimNGrBmwsObJD6lJgIzoDBw7E4cOHsXHjRsTHx8vbhg0bcPToUQwZMuRF6yYiMkzeVYE+fwK2LsDNUCBkJtQqqIQr5oppK0tz7LkQi6ELQ5GRzZEdMk7PPKJjb2+PrVu3okGDBg89vnfvXjRv3hwpKSlQE47oENEziToJhC0Amv8AaNQ9CnLgn9voHxyC9Cwt6pcpium9asKBu7HI1Ed0ihYt+tjpKfGYi4vLs1dKRKQmnlWAlj/fDznZGUD8dahRvdJumNNH7L4yx/5Ld+SurLsp7I1FxuWZg85nn30mz9KJioq695i4/uCDD/D555/nd31ERIZLDIivGwHMbAzcOAo1qlfGTTYCFYcKHr8ej87TDyIqQZ19vojyZeqqWrVquHTpEjIyMuDn5ycfi4iIgLW1NcqWLfvQc8PCwmDoOHVFRM8tLR6Y1waIOqHrft5hOlCpPdToYnSSbAAanZgBHxdbechgCTd7pcsieuHf3888Gdu+vTq/iYmI8p3YddVvM7BqAHBhC7CiD3D3a6D+24CZGdSkrIcjVr5VD71mH8bVO6noNO0gFgyohYpe/AcgqdsLHRhoDDiiQ0QvTJsDbP0EODxNd796H6DVr4C5pdKVPbOYpHT0nn0E56KS4GRjgbn9glDD31XpsogK/8BAIiL6P7EwucVPQPOfADMNEDYPWD0IauTuaINlQ+qipr8LEtOz5QLlXedjlC6LqOBHdPTtTH758mWoCUd0iChfnd8CrBkCdFkIlGwItUrLzMFbC0Ox+0IsLM3NML5zVbQJ9Fa6LKKCW6Nz9epV+Pv7o3v37nB3d9f304iITEv55sA7Jx/uh5WVBljaQk1ED6yZvWti9PJwbDhxC6OWHpMNQbvX1m1CIVILvYPOsmXLMGfOHIwfPx4tWrRA//790bJlS2g0nP0iInrIgyEn5pyudUTLcUBAW6iJlYUGE7tWg7OtJRYdjsAna04iPi0Tw14po3RpRHrTO6W8+eab2Lx5s9xaXqNGDbz77rvw9fXFxx9/jIsXL+r/jkREpuTwVCA5CljeG9g/SXf2joqYa8wwtn1lDHultLz/85bz+GHzWZj4PhZSkWcejilevDg+/fRTGW4WL14s+15VqFABd+/ehZqwezkRFYqWvwJBA8XpgsC2z4EN7wI52VATMzMzfNi8Aj5pWUHen777MsasPokcLcMOGen28vT0dKxcuVJOZR06dAht27bFvHnz5KGBasPFyERU4MSP2UNTdVvQReAp/RrwZvDDU1wqsSwkQoYckXFaVvHEb12qwtpC3T2/SJ0KZHu5GL0ZPHgwPD095Vqdjh074ubNm1i6dKkqQw4RUaEQhwfWHQZ0XQRY2gH/7ADmNAeSoqE2XYL8MKV7dViZa7DpZBQGzjuKlAx1jVCRadE76FSqVAmtW7eGra0tdu/eLds7jBgxgo08iYj0VaEV0Hcj4OAB2LkCtur8+dmiihdm960pm4HuvXhbto6IT2UzUFL51JXYXWVvbw8LCws5X/skcXFxUBNOXRFRoRPdzq3sdWFHxcIi7qLf3BAkpGWhvIejbBnh7mSjdFlkIhLz+xyduXPn5ldtRESmrYjvw/e3fgo4FQfqDFVVj6zqfi5YPqSu7I91PjpJ9scSzUD9itopXRrRPex1xREdIlLSlT26DuhC0CCg+Y+A+TP3W1ZUxJ1UOX0VEZcKd0drzB9QCxU8+fOUVLQY2cSzEBFRwSnREGj6re46ZCawtBuQkQQ1ESM4K9+qiwqejohJykCX6Ydw9lai0mUR6R90xEJksbMqM/Ppi83E2TpDhw7Fjz/+qM/LEhGRmKqqPwroPB+wsAEu/gXMaQEk3ISaiLU5ywbXRVXfInLNjpjO+ic2WemyiPSbutqxYwc++ugj2bCzadOmqFmzJry9vWFjYyMPCjxz5gz27duH06dPy51Yn3zyiRxOUgNOXRGRwbhxFFjSFUiJBRy9gO7LAK9AqIkIOd1mHMKZW4nwcrbBirfqwseFa3ZIud/fz7RGR4QZ0fNq7969uHbtGtLS0uDm5oZq1aqhWbNm6NGjh+q2mzPoEJFBuXsNWNwZiD2vG+VRWX8s4U5yBjpPP4h/YlPgX9QOK4bU5W4sUkfQMUYMOkRkcNLigcs7gUodoFZRCel4c/oBXI9LQzkPBywdXBeu9lZKl0VGpEBORiYiokJgW+ThkCPO3fl7LKDNgVp4Ottg0YA68HCyxoXoZPSZcwSJ6VlKl0UmiEGHiMiQiQagS7oBe8YBy3oB2Zmq2o21aGBtOZJz8mYCBgSHIC1TPWGNjAODDhGRIRNn6jR6T7cj6/xGYO1QQKuFWpRxd8T8/rXgaGOBkKt3MXjBUWRkM+xQ4WHQISIydGIaSzQE1VgCp1YCWz7SdURXicrFnRHcL+heb6xRS44hO0c9YY3UjUGHiEgNyjQBOkwTe0iAIzOA3T9BTWr4u2Jm75qy6/nW09H4YOUJaLXqCWtkYkEnKysL169fx/nz51XXxDPPlClTEBAQgKCgIKVLISLST5VOQMtxuutdPwBh86Em9cu4YUqP6jDXmGHNsZv4fN0pnrxPBU7v7eVJSUlYuHChPCH5yJEj8pRk8amik7mPjw9ef/11DB48WHXBgdvLiUh1dv2km8LqtQZw9oHarAu/iXeWhcvZtyGNSuHjFhXk7xIixbaXjx8/HiVKlJAdzJs0aYK1a9ciPDwcFy5cwMGDB/Hll18iOztbhp3mzZvLVhBERFRAXv4QGPS3KkOO0K5qcfzQoYq8nr7nMn7/+5LSJZGpj+h069YNn332mex59TTp6ekIDg6GlZUV+vfvDzXgiA4Rqd7pNbqWEX51oCaz913BtxvOyOvPWwdgQIOSSpdEKsKTkfXEoENEqnZhK7C4C2DjBPTbDHg8/R+khmbi9ov4bfsFef1jxyroWstP6ZLI1E9G3rlz51MX+BIRUSEq0RDwrQ2kJwALOgJxV6Amo14rg8GNSsnrMWtO4s/jkUqXREbmmYNOx44dERoa+sjjEydOxJgxY/KrLiIi0oeVHdB9KeBeCUiOAhZ0AJKioRZiEfKYFhXQo7afXJw8elk4tp9RT/1khEFn3LhxaNGiBc6dO3fvsV9//RVffPEFNm7cmN/1ERHRf7F1AXqtBor4A3evAAvf0DUGVVHY+bZdZXSoVhzZ2lwMWxyG/ZduK10WmWrQGThwIN5//325++rq1av46aef8M0332DTpk1o2LBhwVRJRERP5+gJ9F4L2LsD0SeBJV2BrHSohUZjhnGdXkKzSh7IzNZi4LyjCL2mznPayAgODPzwww/Ro0cP1KxZEz/++CO2bt2K+vXr5391RESkP9dSupEda2fAuxpgbgU1sTDXYFK3amhY1g1pWTnoOzcEp24mKF0WqZxeu64mTZr02Md/+eUXNGrUCLVq1br32KhRo6Am3HVFREYn4QbgVFzMCUGNRIfz3nMOyyagovP58iF1ZHNQogLbXl6yZEm951kvX74MNWHQISKjlp0BnFoFBHZTVfBJTM9Cj5mHcfJmAjydbLDirbrwdbVTuiwyIDxHR08MOkRktLRaYGFH4PJO4NXPgEYfQE3iUjLRZfpBXIxJhq+rLVYMqQdPZxulyyJjP0eHiIhUQqMByjXTXf89FgiZDTUR01aLBtaGf1E7XI9LQ8/Zh2X4IXoWegUdseA4NTVVrxc8fPgwt5kTERmKOkPvj+RsfA84tRpq4u5kg4UDasPL2QaXYpLl2h0xrUWUr0HnzJkz8Pf3x7Bhw7B582bExsbe+5ho5nnixAn88ccfqFevHrp06QJHRy4aIyIyGI0/BWqK/oO5wOrBwKUdUBOxNmfhwNooam+FUzcT0X9uCFIzs5Uui4wp6MyfPx/bt29HVlYWunfvDk9PT9m4UwQaa2trVKtWDXPmzEHv3r3lQYJiJxYRERkIsQi55S9ApQ6ANgtY1gu4cRRqUrqYAxYMqA0nGwscvXYXQxaEIj0rR+mySAWeeTGyVqvF8ePHERERgbS0NLi5uaFq1aryTzXiYmQiMhnZmcCSLsCNUKD7MsC/LtQmLOIues46jNTMHDQN8MAfParD0pzLTU1RIndd6YdBh4hMSkYykHAdcK8ItTpw6Tb6BofIE5TbV/XG+M5V5cnKZFoS83vXVU5Ojmz3IE5ADgoKwscffyxHdIiISEWsHR4OOdFngOT76y7VoF4ZN0ztUR0WGjOsDY/EZ+tOwcT/zU5PoXfQ+f777/HJJ5/AwcEBxYsXl93Khw8fru+nExGRobl+BJjTXHfWTnqi0tU8k9cqeuC3LlUhBnIWH47AD5vPMezQiwUdsSBZ7KwSfa3Wrl2L9evXY9GiRXLNDhERqZBdUcDcEog6ASzuopvWUpE2gd74seNL8nrGnsv4/e9LSpdEag46YvFxy5Yt790X3ctFy4fIyEio0ZQpUxAQECCn4YiITFLR0kDPVYC1ExBxAFjUCchIgpp0DvLFF60D5PX4bRcwe98VpUsitQYdcV6Ojc3DR29bWlrKLedqJKbdxPlAISEhSpdCRKQc76pArzW6jucRB4GFb6huGqt/g5J4r2k5ef3thjNYFhKhdElkQCz0faKY++zbt688NydPeno63nrrLdjb2997bPVqdZ26SURk8nxqAr3XAgvaA9cP69bs9FwN2KhnJ+qIV8sgOSMb0/dcxserT8LWygJtA72VLovUFHT69OnzyGM9e/bM73qIiEgJxasDvf8E5rcDrBx0a3dURCyl+LhFBRl2Fh2OwOhl4XC0sUDj8u5Kl0YK4zk6PEeHiOi+2POAsy9gZQc10mpz8e7ycKwLj4SNpUb2yapZwlXpsqgAsHs5ERE9u2Ll74cc8e/gwzOA1DiohTg48Jc3A9G4fDGkZ2nRPzgEZ2+pa80R5S8GHSIierx944HNHwDz26oq7IiWEH/0qIGgEi5ITM9Gr9lHcO1OitJlkUIYdIiI6PHKtwTsiwFRJ4F5bYGUO1ALWytzzOoThIpeTridnIGesw8jOjFd6bJIAQw6RET0eKJVRN+NgIMHEC3CThsg5TbUwtnWEvP6B8G/qB2ux6Wh9+wjiE/NVLosKmQMOkRE9PQ1OzLseAIxp4Hg1kByDNTC3dFGLkj2cLLG+egk9AsOQWpmttJlUSFi0CEioqdzK6sLO45eQOxZ3TRWtnpGRnxd7TC/f205wnMsIh5DFoQiIztH6bKokDDoEBHRf3Mrows7Yut53eGAhRXUpLynI+b2C4KdlTn2XryN0cuOI0dr0qermAwGHSIi0r831vDDQPVeUKPqfi6Y3qsGLM3NsPHkLXy29hQ7npsABh0iItKf1f2WP0iOBZZ0AxJuQi0ali2GCV2qwcwMWHIkAuO2nle6JCpgDDpERPR8/hwJnN8EBLcE4q9DLVq95IXvO1SR13/s+gcz91xWuiQqQAw6RET0fFqOA1xKAHevAsGtgHj1dA3vVssPHzYvL6+/23QWy4+qJ6jRs2HQISKi51PEV7dA2aUkEH8NmNtKF3pUYujLpTG4USl5/fGqE9hyKkrpkqgAMOgQEdHzc/YB+m0CXEsDCRG6c3birkAtHc/HtKiAzjV9IDZgjVpyDAcuqedARNIPgw4REb0YJ2/dyE7RskDCdWDdcF1DUJWEHbFep1klD2TmaDF4QSjORLIJqDFh0CEiohfn5AX03QCUaQp0mCYSBNTCwlyDiV2roXZJVyRnZKPv3CO4HpeqdFmUTxh0iIgofzh6Aj1XAkX87j+mkpEdG0tzzOhdE+U9HBGTlIE+c4/gbop6Tn+mJ2PQISKignF2PTC/HZCpjtER0SIiuH8QvJxtcDk2BQPnH0V6FltFqB2DDhER5b/0RODPUcCV3cCKvkBOFtTAy9kW8/rXgpONBUKv3cXIJceQnaNVuix6AQw6RESU/2ycgK6LAAsb4OJW3eGCWnUEhnIejpjVJwhWFhpsOxONL/48zVYRKsagQ0REBcO/HvDmPMDMHDi+BNj2uWrW7NQq6YqJXarKNdWLD0dg8t+XlC6JnhODDhERFZzyzYF2U3TXBycD+36DWrSo4oWv2lSS179uu8DTk1WKQYeIiApW1W7A69/prnd8DVzeBbXoU68Ehr5SWl6PWX0SO8/FKF0SPSMGHSIiKnj1RgAN3gVqDQFKNIKafNisPDpWK44cbS6GLQpD+PV4pUuiZ2CWa+IrrBITE+Hs7IyEhAQ4OTkpXQ4RkfHK+3WjosME82TlaNE/OAR7L96Gq70VVg+thxJu9kqXZdIS9fz9zREdIiIqHCLg5IWc7Exg3Qgg8hjUwNJcg6k9a6BycSfEpWSi95wjiE3KULos0gODDhERFb7dPwHHFgALOwG31bGjycHaAnP6BsHX1RYRcalyhCclI1vpsug/mGzQmTJlCgICAhAUFKR0KUREpqf+24BXIJB6G1jQAUiMhBq4O9pgXr9acvrq5M0EDF0UJqe1yHBxjQ7X6BARKSM5FpjTDIj7ByhWEei3CbBzhRoci7iL7jMPIy0rBx2rF8evbwbKTuhUeLhGh4iIDJtDMaDXGsDRC4g9CyzuAmSmQA2q+blgSo9qMNeYYXXYTYzbel7pkugJGHSIiEg5Lv5Az9WATRHgxhFg9WCoxasVPPB9h8ry+o9d/2DegatKl0SPwaBDRETK8ggAui8HHDyA2m9BTboE+eHdJuXk9VfrT2PzyVtKl0T/wqBDRETK86sNvH0cKNkQajPqtTLoVstPHhP09rJwHLkSp3RJ9AAGHSIiMgyWtvevY84Ch6dDDcQi5G/bVUKTih7IzNZi4LwQXIhOUros+j8GHSIiMizJMcCc5sDmD4HQeVADC3MNfu9WDdX9iiAxPRt95hzBrYQ0pcsiBh0iIjI4Du5A0ADd9YZ3gPOboQa2VuaY3ScIpYvZ41ZCOvrOCUFCWpbSZZk8Bh0iIjI8r34OVO0J5GqBFX2BiMNQAxd7K8zrXwvujtY4H52EQfOPIj0rR+myTBqDDhERGR5x+F6biUDZZkB2OrC4MxBzDmrg42KH4H614GhtIRcmj14eLjufkzIYdIiIyDCZWwBvzgWK1wTS44GFHYGEm1CDAG8nTO9VA5bmZth0MgrfbjgDE29EoBgGHSIiMlxW9rozdoqWBVxLAdYOUIt6Zdzwa+eq8jr4wFVM231Z6ZJMkoXSBRARET2VfVGgz3rA1gWwtIGatA30RkxiOsZuPIuftpyDh5M1Olb3Ubosk8IRHSIiMnxOXg+HnHMbAa06FvkObFgKgxqWlNcfrTqBA5duK12SSWHQISIiddn+NbC0O7DxPcjjiFVgTIuKaPWSF7JycjFkYSgu8kDBQsOgQ0RE6uIVKLZlAaFzgd0/Qw00GjP8+mYgavq7ICk9G33nhiAmKV3pskwCgw4REalLpfZAy3G6613fA6HBUAMbS3PM6F0TJd3scTM+DQOCjyI1M1vpsowegw4REalPrUFAw/d11xveBc5tghq42lthbt8guNhZ4uTNBIxacoxn7BQwBh0iIlKnVz8Dqv3/9OSV/YCIQ1CDEm72mNWnJqwsNNh+NgbfrD/NM3YKEIMOERGp9/Tk1hOBcs2B7Awg5gzUooa/KyZ00Z2xM+/gNczed0XpkowWgw4REan79OROc4EeK4Ga/aEmLat44ZOWFeT1d5vOYsupW0qXZJQYdIiISN2s7ICyTe7fT4sH0u5CDQY1LIWedfzkLvm3l4bjWIQ66lYTBh0iIjIeiZHA3BbA0h666SwDZ2Zmhq/aVMKrFdyRka3FwHlHce1OitJlGRUGHSIiMh5iNCfhBnBtP/DnKFUcKGhhrsHv3aqhkrcT7qRkot/cEMSnZipdltFg0CEiIuPhEQC8GQyYmQMnlqrmQEF7awvM6RsEb2cbXL6dgsHzQ5GepY4WF4aOQYeIiIxLmdeA1uPvHyh4YjnUwMPJBnP71YKjtQWOXI3DhytPQMszdl4Ygw4RERmfGn2BeqN01+uGA9cOQA3Kezpias8asNCY4c/jkRi/7YLSJakegw4RERmnJl8DFdsCOZnApg8BrRZq0KCsG77vWEVeT955CctCIpQuSdUslC6AiIioQGg0QIfpgI0z0PgT3X2V6FzTFzfiUjHp70v4ZM0peDnbolG5YkqXpUrq+VsnIiJ6njN22k0GnLyhNu82LYcO1YrLXljDFoXhXFSi0iWpEoMOERGZjlOrgHUjVLHtXJyx8+MbVVC7pCuSM7LltvPoxHSly1IdBh0iIjIN8RHA6iHAsQXArh+hBtYW5pjRqyZKFbPHrYR0GXZE6CH9MegQEZFpKOIHtPpFd737R+D4MqiBs50lgvvWQlF7K5y5lYiRi8OQnaOOhdWGgEGHiIhMa9t5/Xfubzu/uh9q4FfUDrP61IS1hQY7z8fiyz9PI1cF02+GgEGHiIhMy2tfAgHtAG0WsLQ7cPsi1KCanwsmdq0KMzNg0eEIzNx7WemSVIFBh4iITHPbefGaQHo8sOhN1XQ7b17ZC5+2rCivv990DhtP3FK6JIPHoENERKbH0hbotkS3biegLWDtDLUY0KAk+tT1l9fvLg9H6DV1hDSlMOgQEZFpcnAHhuwFmn6jqsMExbbzL9pUQpOK7sjM1mLw/KO4HpeqdFkGSz1/s0RERPnNtsj966x04MyfUANzjRkmdq2GAC8n3EnJRP/gECSkZSldlkFi0CEiIsrOAOa3A5b3Uk23c3trC8zuWxMeTta4GJOMEYvDkMVt549g0CEiIrKwBvxqP9Dt/CDUQPTAmt0nCLaW5th78Ta+WMdt5//GoENERCS89hVQobWu27nYdh6nju3blYs7Y1K3anLb+ZIjEZi194rSJRkUBh0iIiJBLEjuOAPwqgqkxQGLu6hm23nTAI/72843n8XW01FKl2QwGHSIiIjyWNkD3ZYCTsWB2xeA5b2BnCzVbDvvWcdP9it9e+kxnLyRoHRJBoFBh4iI6EFOXkD3ZYCVA3DrBHDnH6hl2/lXbSqhYVk3pGdpMWBeCCLj02DqVB904uPjUbNmTVStWhWVK1fGzJkzlS6JiIjUzrMK0GUBMHAH4F4BamFhrsGUHtVRzsMBMUkZctu5qXc7N8tV+fLsnJwcZGRkwM7ODikpKTLsHD16FEWLFtXr8xMTE+Hs7IyEhAQ4OTkVeL1ERKRS4pwdSxuowY27qWg/ZT9uJ2eicflimNm7pgxBxkTf39+q/6rNzc1lyBFE4BG5TeXZjYiIDM0/O4FJVYHIY1ADHxc7GW6s/9/tfOzGszBVigedPXv2oE2bNvD29pbzi2vXrn3kOVOmTEGJEiVgY2OD2rVr48iRI49MXwUGBsLHxwcffPAB3NzcCvErICIio3d4GpB0C1jcFUi4CbV0O/+tS1V5HXzgKhYcvApTpHjQEdNNIqSIMPM4y5Ytw+jRo/Hll18iLCxMPrdZs2aIiYm595wiRYrg+PHjuHLlChYvXozo6OhC/AqIiMjodZwJuAcAyVHAki5ARjLUoGUVL3zQrLy8/mr9Gey+EAtTo3jQadGiBcaOHYsOHTo89uPjx4/HoEGD0K9fPwQEBGDatGlyqmrOnDmPPNfDw0MGob179z7x/cT0lpjXe/BGRET0VDZOup1Y9sWAqJPAqoGANgdqMOyV0uhYvThytLkYsSgMF6KTYEoUDzpPk5mZidDQUDRp0uTeYxqNRt4/eFB3PLcYvUlK0v2liQVJYiqsfHlden2cH374QS5eyrv5+voWwldCRESqV8QP6LoEsLABLmwG/vocamBmZoYfOlZBrRKuSMrIljuxbidnwFQYdNC5ffu23FUlRmoeJO5HRelOfbx27RoaNmwoR3LEnyNHjkSVKlWe+JpjxoyRgSjvdv369QL/OoiIyEj4BgHtp+quD00Bzq6HGlhbmGNarxrwL2qHG3fTMGRBKNKz1DEi9aIsoHK1atVCeHi43s+3traWNyIioudSuaOuD5a4lW0GtXC1t5INQDv+sR+h1+7io1UnMKFLVTniY8wMekRH7J4S28f/vbhY3Pf09FSsLiIiMnEN3wPaTQEsrKAmZdwdMLVnDZhrzLAuPBKTdlyCsTPooGNlZYUaNWpgx44d9x7TarXyft26dRWtjYiITJgYBckbCRGLknf9CCSpo5Fm/TJu+LZdZXn92/YL+PN4JIyZ4lNXycnJuHTpfqIUW8TFVJSrqyv8/Pzk1vI+ffrINg9immrChAlyS7rYhUVERKS4rZ/oztm5sAXou1HXGNTAda/th8uxyZi17wreX3EcPi62qO7nAmOkeAuIXbt2oXHjxo88LsJNcHCwvJ48eTLGjRsnFyCLnlaTJk2SBwfmB7aAICKiFyLW6sxqAqTeASq0BjrPBzTmMHQ52lwMWXAU28/GwM3BCmuH15cnKquFvr+/FQ86SmPQISKiFxZxCJjXBsjJBOqOAJp9BzVIychGp2kHcfZWIip4OmLl0HpwsFZ8skcvJtPr6nmJk5jFAYRBQUFKl0JERGrnV+f+tvODk4GQ2VADe2sLzOpTE24O1jgXlYS3lxyTIz3GhCM6HNEhIqL8snscsHMsYGYOdF8OlL1/4K0hOxZxF11mHEJmthaDG5XCJy0rwtBxRIeIiKiwNXofCOwOmFsBOeo5fbianwt+eTNQXs/YcxnLQiJgLBh0iIiI8ovYct5mIjDob6BCK6hJ20BvvP1aWXn96ZpTOHT5DowBgw4REVF+EocIegTcv59wA8hMgRq806QsWr/khWxtLt5aGIqrt9VR99Mw6BARERWUm6HAjMbA6sGq6HZuZmYmp7ACfZwRn5qFAfNCkJCWBTVj0CEiIiooOVlAejxwbgOw7QuogY2lOWb2rgkvZxv8E5uCEYvDkJ2jhVox6BARERXWtvNji6AG7k42MuzYWppj78Xb+GbDGagVgw4REVFBqtIJePlj3fWGd3SHC6pA5eLOmNBVdDcH5h+8hgUHr0KNTDbo8MBAIiIqNC9/BFRsqzs5eVlPIP461KBZJU982KyCvP5q/Rnsu3gbasMDA3lgIBERFQax82pOMyDqJFClM/DGTKhBbm4u3lt+HKuP3YSTjYXsiVWqmIPSZfHAQCIiIoMiupp3XQJU6wW0/g1qYWZmhu87VkE1vyJITM/GwHlHkZCqnp1YDDpERESFpYgv0G4yYK38iMiz7sSa0asmvJ1tcPl2CoaraCcWgw4REZESxMqRvb8Cp9dADYo5WmNmH91OrH2XbuNblezEYtAhIiJSwskVwI5vgDVDgchwqEElb2f81qWqvJ538BoWHroGQ8egQ0REpIRKHYEyTYDsNGBpdyApGmrQvLInPmhWXl5/+edpHLhk2DuxGHSIiIiUYG4BvDEbKFoWSLyp23aerY6O58NeKY32Vb2Ro83F0EVhuGLAPbEYdIiIiJRiWwTovgywcQZuHAHWv6Nbu6OCnVg/vvESqvoWkb2wDLknFoMOERGRkoqWBt4MBszMgeOLda0i1LITq3cN2RPrcmwKRi45ZpA7sUw26PBkZCIiMhilXwWafQ9oLAAr9Ww9d3e83xNrz4VY/LD5HAwNT0bmychERGQIxK/jO5cAt7JQm00nb2HYojB5/dMbVdAlyK/A35MnIxMREamJ6J75YMhJjQNSDHtHU56WVbzwbpNy8vqztadw5EocDAWDDhERkaGJvQDMfFW37TwrHWow6rUyaPWSF7JycvHWwlBcj0uFIWDQISIiMkRpccD1w8D6UarZifVLp0BULu6EuJRMDJp/FMkZ2UqXxaBDRERkcIqVA96cp9uJdWIZsPcXqIGtlblcnCzaRZyLSsI7S8Oh1Sob0hh0iIiIDFHpxkDLcbrrv8cCp9dCDbycbTGjVw1YWWiw/Ww0fvnrvKL1MOgQEREZqqABQO2huus1bwE3dTubDF01Pxf8/MZL8vqPXf/g73PKtbewUOydiYiI6L81+0637fzSNmDrp0C/TbodWgaufbXiuBCdhNikDNQv46ZYHTxHh+foEBGRoUtPBP76DGjyFWDnCrXQanNlJhMLlZX6/c0RHSIiIkNn4wS0nQS10WiUH3niGh0iIiK1CQ0G/v5O6SpUwcKUe12JW05OjtKlEBER6S8yHFj/tu66aBkgsIvSFRk0rtHhGh0iIlKbHd8Ae38FzK2AvhsB31owNYnsdUVERGSkGn8GVGgN5GTq2kTERyhdkcFi0CEiIlIbjQboOAPwrAKkxAJLugEZyUpXZZAYdIiIiNTIyh7othSwdweiTwGrB4n93EpXZXAYdIiIiNTK2QfouhiwsAWK11DFQYKFzWR3XRERERkF3yBg1DHAyUvpSgwSR3SIiIjU7sGQI9bqRJ9WshqDwqBDRERkLJKigLnNgUVvAtkZSldjEBh0iIiIjIW1I2BfDHhjNmBhrXQ1BoFrdIiIiIxpJ1avNUpXYVA4okNERERGi0GHiIiIjBaDDhERERktkw06onN5QEAAgoKClC6FiIiICgi7l7N7ORERkeqwezkRERGZPAYdIiIiMloMOkRERGS0GHSIiIjIaDHoEBERkdFi0CEiIiKjxaBDRERERotBh4iIiIwWgw4REREZLQuYuLyDocUJi0RERKQOeb+3/6vBg8kHnaSkJPmnr6+v0qUQERHRc/weF60gnsTke11ptVpERkbC0dERZmZmitQgGouGhITAWBn616dUfYX1vgXxPvn1mvnxOs/zGuJfguIfN9evX2ePO4UY+s8FY//6ghSsL7/eW8QXEXK8vb2h0Tx5JY7Jj+iI/3N8fHwUrcHc3Nyof9ga+tenVH2F9b4F8T759Zr58Tov8hri8wz5v01jZug/F4z96zNXsL78fO+njeTk4WJkAzB8+HAYM0P/+pSqr7DetyDeJ79eMz9ex9D/+yLT/Hsz9K9vuIL1FfZ7m/zUFRGZHjF1Jf4lmJCQYND/6iaiF8cRHSIyOdbW1vjyyy/ln0Rk3DiiQ0REREaLIzpERERktBh0iIiIyGgx6BAREZHRYtAhIiIio8WgQ0REREaLQYeI6D906NABLi4u6NSpk9KlENEzYtAhIvoPb7/9NubPn690GUT0HBh0iIj+wyuvvCIb/xKR+jDoEJGq7dmzB23atJEdjM3MzLB27dpHnjNlyhSUKFECNjY2qF27No4cOaJIrURU+Bh0iEjVUlJSEBgYKMPM4yxbtgyjR4+WLR/CwsLkc5s1a4aYmJh7z6latSoqV678yC0yMrIQvxIiKghsAUFERkOM6KxZswbt27e/95gYwQkKCsLkyZPlfa1WC19fX4wcORIff/yx3q+9a9cu+RorV64skNqJqGBwRIeIjFZmZiZCQ0PRpEmTe49pNBp5/+DBg4rWRkSFg0GHiIzW7du3kZOTAw8Pj4ceF/ejoqL0fh0RjN58801s2rQJPj4+DElEKmKhdAFERIZu+/btSpdARM+JIzpEZLTc3Nxgbm6O6Ojohx4X9z09PRWri4gKD4MOERktKysr1KhRAzt27Lj3mFiMLO7XrVtX0dqIqHBw6oqIVC05ORmXLl26d//KlSsIDw+Hq6sr/Pz85NbyPn36oGbNmqhVqxYmTJggt6T369dP0bqJqHBwezkRqZrY9t24ceNHHhfhJjg4WF6LbeHjxo2TC5DFmTmTJk2S286JyPgx6BAREZHR4hodIiIiMloMOkRERGS0GHSIiIjIaDHoEBERkdFi0CEiIiKjxaBDRERERotBh4iIiIwWgw4REREZLQYdIipQr7zyCt55551n/rzPP/8cgwcPhtpMmzYNbdq0UboMIvo/Bh0iMjiiVcPEiRPx6aef3nusb9++MDMzw1tvvfXI84cPHy4/Jp5TkNLT0+V7VKlSBRYWFmjfvv0jz+nfvz/CwsKwd+/eAq2FiPTDoENEBmfWrFmoV68e/P39H3rc19cXS5cuRVpa2kPhY/HixbKBZ0HLycmBra0tRo0ahSZNmjyxY3r37t1lPy0iUh6DDhEVqo0bN8LZ2RmLFi164nNEmHnc9E/16tVl2Fm9evW9x8S1CDnVqlV7ZMpsxIgR8ibez83NTU6HPdjeLyMjAx999JF8TWtra5QpUwazZ89+Yl329vaYOnUqBg0aBE9Pzyc+T9T+559/PhTIiEgZDDpEVGjEyEu3bt1kyOnRo8djnxMXF4czZ86gZs2aj/24mBqaO3fuvftz5sxBv379HvvcefPmySmmI0eOyKmw8ePHy9GiPL1798aSJUvk6MvZs2cxffp0ODg4vPDXKWrPzs7G4cOHX/i1iOjFWLzg5xMR6WXKlClyzc369evx8ssvP/F5ERERctTF29v7sR/v2bMnxowZg2vXrsn7+/fvlyNAu3bteuS5YqTmt99+k+t3ypcvj5MnT8r7YkTmwoULWL58ObZt23ZvGqpUqVL58rXa2dnJUaS8GolIOQw6RFTgVq5ciZiYGBlKgoKCnvrcvOkeGxubx368WLFiaNWqFYKDg2UgEtdiWupx6tSpI0NOnrp16+LXX3+Va23Cw8Nhbm7+xNBVqVKle0GlYcOG2Lx5M56FWMuTmpr6TJ9DRPmPQYeICpxYPyN2IolpJjGt82D4+Le80HL37l0Zap40fSXW3uSNFD0PEUSeZtOmTcjKytLruU+agntS/URUeLhGh4gKXOnSpbFz506sW7cOI0eO/M/nOjk5yXU6T9K8eXNkZmbKINKsWbMnPu/fa2QOHTqEsmXLypEcsUVcq9Vi9+7dj/1cseNLLE4Wt+LFi+NZ/PPPP3I32L8XSBNR4WPQIaJCUa5cORl2Vq1a9dQDBDUajVwzs2/fvic+RwQVsXhYhCFx/bT1PqNHj8b58+flouPff/8db7/9tvxYiRIl0KdPHzk6tHbtWly5ckWu8xHrdp5GvKeY9hIjNgkJCfJa3B4kztAR631EaCMiZXHqiogKjVgQ/Pfff8ut3yKgiPUyjzNw4EC5YPjnn3+WwedxxKjPfxG7qsSan1q1asn3EyHnwdOWxVbxTz75BMOGDcOdO3fkNnVx/2latmz50CLjvFGbB7eti1Al6ici5ZnlPvjdSURkAMSPpdq1a+Pdd9+V29GfhwhTVatWxYQJE1CYTp8+jVdffVXu6hI7r4hIWZy6IiKDIxYrz5gxQ55Foza3bt3C/PnzGXKIDASnrojIIInRGHFTmye1hiAiZXDqioiIiIwWp66IiIjIaDHoEBERkdFi0CEiIiKjxaBDRERERotBh4iIiIwWgw4REREZLQYdIiIiMloMOkRERGS0GHSIiIgIxup/ajjMbSXw+dIAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(keff, power_model_1d_1,label='true box')\n", "plt.plot(keff2, power_model_1d_2,label='cubic box',ls='--')\n", "\n", "plt.xscale(\"log\")\n", "plt.yscale(\"log\")\n", "plt.xlabel(\"k (Mpc-1)\")\n", "plt.ylabel(\"P(k) (Mpc3)\")\n", "plt.legend()" ] }, { "cell_type": "markdown", "id": "ab32078f", "metadata": {}, "source": [ "You can see that, even just for the matter power without any FoG, beam, etc, the effect of assuming the wrong k-sampling in your forecast will already be massive." ] }, { "cell_type": "code", "execution_count": null, "id": "d87e4cef", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "meer21cm", "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.10.16" } }, "nbformat": 4, "nbformat_minor": 5 }