{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Computing power spectra\n\n[![Binder](https://mybinder.org/badge_logo.svg)](https://binder.flatironinstitute.org/v2/user/fvillaescusa/Quijote?filepath=/Tutorials/Power_spectrum.ipynb)"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport readgadget\nimport readfof\nimport MAS_library as MASL\nimport Pk_library as PKL\nimport redshift_space_library as RSL","execution_count":1,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We start by computing the matter power spectrum of a snapshot"},{"metadata":{"trusted":true},"cell_type":"code","source":"snapshot = '/home/jovyan/Data/Snapshots/s8_p/397/snapdir_004/snap_004' #location of the snapshot\n\n# density field parameters\ngrid    = 512    #the density field will have grid^3 voxels\nMAS     = 'CIC'  #Mass-assignment scheme:'NGP', 'CIC', 'TSC', 'PCS'\nverbose = True   #whether to print information about the progress\n\n# power spectrum parameters\naxis = 0    #axis along which redshift-space distortions have been placed. In real-space this parameter doesnt matter\nthreads = 1 #number of openmp threads to compute the power spectrum","execution_count":2,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"First, lets read the particle positions:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# read the redshift, boxsize, cosmology...etc in the header\nheader   = readgadget.header(snapshot)\nBoxSize  = header.boxsize/1e3  #Mpc/h\nNall     = header.nall         #Total number of particles\nMasses   = header.massarr*1e10 #Masses of the particles in Msun/h\nOmega_m  = header.omega_m      #value of Omega_m\nOmega_l  = header.omega_l      #value of Omega_l\nh        = header.hubble       #value of h\nredshift = header.redshift     #redshift of the snapshot\nHubble   = 100.0*np.sqrt(Omega_m*(1.0+redshift)**3+Omega_l)#Value of H(z) in km/s/(Mpc/h)\n\n# read positions of the dark matter particles\npos = readgadget.read_block(snapshot, \"POS \", [1])/1e3 #positions in Mpc/h","execution_count":3,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Second, lets compute the density field:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# define the matrix hosting the density field\ndelta = np.zeros((grid,grid,grid), dtype=np.float32)\n\n# construct 3D density field\nMASL.MA(pos, delta, BoxSize, MAS, verbose=verbose)\n\n# compute the overdensity field\ndelta /= np.mean(delta, dtype=np.float64);  delta -= 1.0\n\n# print some information\nprint('%.3f < delta < %.3f'%(np.min(delta), np.max(delta)))\nprint('< delta > = %.3f'%np.mean(delta))","execution_count":4,"outputs":[{"output_type":"stream","text":"\nUsing CIC mass assignment scheme\nTime taken = 5.516 seconds\n\n-1.000 < delta < 1044.773\n< delta > = -0.000\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"Third, compute the power spectrum"},{"metadata":{"trusted":true},"cell_type":"code","source":"# compute power spectrum\nPk = PKL.Pk(delta, BoxSize, axis, MAS, threads, verbose)\n\n# Pk is a python class containing the 1D, 2D and 3D power spectra, that can be retrieved as\n\n# 1D P(k)\nk1D      = Pk.k1D\nPk1D     = Pk.Pk1D\nNmodes1D = Pk.Nmodes1D\n\n# 2D P(k)\nkpar     = Pk.kpar\nkper     = Pk.kper\nPk2D     = Pk.Pk2D\nNmodes2D = Pk.Nmodes2D\n\n# 3D P(k)\nk       = Pk.k3D\nPk0     = Pk.Pk[:,0] #monopole\nPk2     = Pk.Pk[:,1] #quadrupole\nPk4     = Pk.Pk[:,2] #hexadecapole\nPkphase = Pk.Pkphase #power spectrum of the phases\nNmodes  = Pk.Nmodes3D","execution_count":5,"outputs":[{"output_type":"stream","text":"\nComputing power spectrum of the field...\nTime to complete loop = 7.10\nTime taken = 12.18 seconds\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"Lets see how the 3D matter power spectrum looks like:"},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nplt.xscale('log')\nplt.yscale('log')\nplt.xlabel(r'$k~[h{\\rm Mpc}^{-1}]$')\nplt.ylabel(r'$P(k)~[(h^{-1}{\\rm Mpc})^3]$')\nplt.plot(k, Pk0)\nplt.show()","execution_count":6,"outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 640x480 with 1 Axes>","image/png":"iVBORw0KGgoAAAANSUhEUgAAAkIAAAG9CAYAAAD5ixlRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABII0lEQVR4nO3dd3hUZd7G8XvSSUgPBBISQglgKAmEBFBAqlgWBdYuXREUXVZ0Lbsq++7q2rExgoI0CyKKrLuuiCJdOoQWqiCEkkAI6aTNzPsHOoq0lEnOJPP9XFcuZ86ceeaXcCQ352kmm81mEwAAgAtyM7oAAAAAoxCEAACAyyIIAQAAl0UQAgAALosgBAAAXBZBCAAAuCyCEAAAcFkEIQAA4LI8jC7AmVmtVh0/flz+/v4ymUxGlwMAAMrBZrMpLy9PERERcnO7/D0fgtBlHD9+XFFRUUaXAQAAKiEtLU1NmjS57DkEocvw9/eXdO4HGRAQYHA1AACgPHJzcxUVFWX/PX45BKHL+KU7LCAggCAEAEAtU55hLQyWBgAALosgBAAAXBZBCAAAuCyCEAAAcFkEIQAA4LIIQgAAwGURhAAAgMsiCAEAAJdFEAIAAC6LIAQAAFwWQegizGaz4uLilJSUZHQpAACgGplsNpvN6CKcVW5urgIDA5WTk8NeYwAA1BIV+f3NHSEDnC2x6IN1h7XreI7KLFajywEAwGWx+7wBdhzL0TOLdkqSfL3c1aFJoDpFB6tjdLA6RQcptL63wRUCAOAaCEIGMJmkHrFhSjmSrbziMq07mKV1B7PsrzcN9VXHqCB1ahqsjlHBatPYX57u3LwDAMDRGCN0GdU9RshitenAyXxtPXJGW46c0dYj2dp/Mv+C83w83dShSZA6RgepU3SwOkUHq4E/d40AALiYivz+JghdhhGDpXPOliolLfvncJStlCNnlFtUdsF5TYLr/dyddi4cXdU4QF4e3DUCAIAg5CDOMGvMarXpYGa+thzO1ta0M9pyOFv7Tubp939q3h5uah8ZqE5Nz40z6hgdrPAAH0NqBgDASAQhB3GGIHQxuUWl2p6W83N32hltTctWdmHpBedFBtVTgr07LUhxEQHy9nA3oGIAAGoOQchBnDUI/Z7NZtOhzAJtOZJtH2u0Nz1X1t/9yXp5uKldRMDPs9OC1alpkBoH1jOmaAAAqglByEFqSxC6mPziMm1Py9bWtGxtOXzurlFWQckF5zUK8FGnpkH28UZtIwLl48ldIwBA7UUQcpDaHIR+z2az6fDpQvsdoy1HzmhPep4sv7tt5OluUtuIQPsg7JYN68tkkn1Mks0m2WQ7b4zSL8d+fXzu86RfHp97ZG/jl/N+d47JJF3VKECBvp7V9WMAALgAgpCD1KUgdDGFJWXafjTHHo62HjmjzPwL7xrVJD8vd426ppnu69FMQb5ehtYCAKidCEIOUteD0O/ZbDalZZ39eXbauen7x7PPymSSJNPP/5VM0m8enzv+81OZfn7B/rrp13N+fa/J/li/ee/ZEouO5xRJkvy9PTS6ezON7t5MgfW4QwQAKD+CkIO4WhAyms1m0ze7MvTGd/u0Jz1PkhTg46H7ezbXyGuaqb43C6EDAK6MIOQgBCFjWK02fb0zXa9/t08Hfl5pO9jXU/f3bKERVzeVrxeBCABwaQQhByEIGctitem/24/rze/262BmgSQp1M9LD/RqoXu6NFU9L2a3AQAuRBByEIKQcyizWPXvlON6c+l+HckqlCQ18PfWg71a6K7kaKb7AwDOQxByEIKQcym1WLVwy1G9tfSAjmWflXRuHaTxfVrq9s5NWDUbACCJIOQwBCHnVFJm1YLNaZry/QGd+HmWWWRQPT3Up6VuTWwiT3c2nwUAV0YQchCCkHMrLrPokw1pMi87oJN5xZKkqJB6+lOfWA3uGCkPAhEAuCSCkIMQhGqHolKLPlp/RFOXH7AvCBkT6qsJ/WJ1c3yk3N1MV2gBAFCXEISqyGw2y2w2y2KxaN++fQShWqKwpEwfrD2saSt+1JnCUklSszA/DezQWNe1baS2EQH2xRwBAHUXQchBuCNUO+UXl2nODz/pvZUHlXO21H48Mqie+seF67q24UqOCaHrDADqKIKQgxCEare8olJ9m5qhJbsytGLfKZ0ttdhfC/L1VN8250JRz9gGrEkEAHUIQchBCEJ1x9kSi1YfyNSSXen6bneGvetMknw83dQztoGua9tIfds0VLBfxTd7LbNYdTy7SEeyCpWeW6RerRsorL63I78FAEA5EYQchCBUN5VZrNp0+IyW7MrQN7vS7WsSSZK7m0nJMSG6rm24+seFq0mwr/21/OIyHT5doLSsQh0+XajDWYX2x8eyz8pi/fV/pRYN/LTwgWsU6MuGsQBQ0whCDkIQqvtsNptST+Rqya4MLUnN0O4Tuee9Htc4QF4ebjqSVaisgpLLtuXl4aao4Ho6U1iqrIISdWseqjmjk+XlwVgkAKhJBCEHIQi5niOnC7UkNV1LUjO06acsWX/3f0eIn5eiQnzVNMRX0SG+ig79+XGor8L9feTmZtLuE7m6bdpa5ReXaUinSL12Wzyz1QCgBhGEHIQg5NpO5xdr9YFMebm7KTrUV1EhvgrwKV9X14p9pzR69kZZrDY90q+VJvSLreZqAQC/IAg5CEEIVfHx+iP66xc7JEmTb4/XkE5NDK4IAFxDRX5/M3gBqCZ3d4nWuGtbSJKe+Hy71v542uCKAAC/RxACqtHjA1rrpg6NVWqxaewHm3TgZL7RJQEAfoMgBFQjNzeTXrstXp2ig5RbVKZRszcoM7/Y6LIAAD8jCAHVzMfTXdOHd1Z0iK/Sss7qvjmbVPSbVa4BAMYhCAE1ILS+t2aNSlJgPU+lpGXrkfkpsv5+bj4AoMYRhIAa0qJBfb03LFFe7m76eme6Xly8x+iSAMDlEYSAGtSleaheua2DJOm9lQf14brDBlcEAK6NIATUsFsSIvVo/1aSpGf/vVPL9pw0uCIAcF0EIcAAD/VpqVsTm8hqkx78aItW7880uiQAcEkEIcAAJpNJ/xrcXr1aN9DZUotGz9mopbszjC4LAFwOQQgwiJeHm94dlqgBbcNVUmbV2A8266vtJ4wuCwBcCkEIMJC3h7um3N1JtyREqMxq08PztujzzUeNLgsAXAZBCDCYp7ubJt+eoDuTomS1SY8u2MZsMgCoIQQhwAm4u50bMzTy6hhJ0tOLdmrGqoPGFgUALoAgBDgJNzeTJg2M0wO9zu1Y/9xXu/XW0v2y2ViBGgCqC0EIcCImk0mPD2htX2do8rf79NLivYQhAKgmHkYXAOB8JpNJD/eNVT0vdz331W5NW/GjikotevYPcXJzM9nPKymzKr+4TPlFZcovLlNhSZlaN/KXv4+ngdUDQO1CEAKc1H09msvH011PL9qp2T/8pBX7TslitZ0LP8VlKimzXvCesPpeevnWDurTJtyAigGg9qFrDHBiQ7s21Wu3xcvNJB3KLNCRrEJlFZScF4Lqebqrgb+3Qvy8lJlfotGzN+npRTt0tsRiYOUAUDuYbAw+uKTc3FwFBgYqJydHAQEBRpcDF3b4dIEOny6Un7eH/H085OftofreHvLzcpeH+7l/zxSVWvTKN3v1/upDkqQWDfz05p0d1S4y0MjSAaDGVeT3N0HoMghCqI1W7T+lRz/dppN5xfJ0N2li/9a6v2dzuf9mfBEA1GUV+f1N1xhQx/SIbaBv/txT17dtpFKLTS8t3qO7p6/TseyzRpcGAE7HJYJQYWGhmjZtqscee8zoUoAaEeznpalDO+nlP3aQr5e71h/K0vVvrNS3qWzsCgC/5RJB6Pnnn1fXrl2NLgOoUSaTSbcnRenrCT3UMTpIeUVlevCjzVq1/5TRpQGA06jzQWj//v3as2ePbrjhBqNLAQzRNNRPC8Z2000dGqvUYtPYDzYrJS3b6LIAwCk4dRBauXKlBg4cqIiICJlMJi1atOiCc8xms2JiYuTj46MuXbpow4YN573+2GOP6YUXXqihigHn5OHuptdvT1CP2DAVllg0atYGHTiZZ3RZAGA4pw5CBQUFio+Pl9lsvujr8+fP18SJEzVp0iRt2bJF8fHxGjBggE6ePClJ+ve//61WrVqpVatWNVk24JS8PNw0bWii4qOCdKawVMPe38AAagAur9ZMnzeZTPriiy80aNAg+7EuXbooKSlJU6ZMkSRZrVZFRUXp4Ycf1pNPPqmnnnpKH374odzd3ZWfn6/S0lI9+uijevbZZy/6GcXFxSouLrY/z83NVVRUFNPnUadkFZTotmk/6MdTBWrewE+fjbtaIX5eRpcFAA7jEtPnS0pKtHnzZvXr189+zM3NTf369dPatWslSS+88ILS0tL0008/6dVXX9WYMWMuGYJ+OT8wMND+FRUVVe3fB1DTQvy89MG9XRQR6KODpwo0atYG5ReXGV0WABii1gahzMxMWSwWhYefv6dSeHi40tPTK9XmU089pZycHPtXWlqaI0oFnE5EUD3NvbeLgn09te1ojsZ9sFnFZWzJAcD1uMymqyNHjrziOd7e3vL29q7+YgAn0LJhfc0alay7p6/T6gOZ6vPqCrWPDFSrRv5q08hfrRv5q2mIr30LDwCoi2ptEAoLC5O7u7syMs5fIC4jI0ONGjUyqCqgdkmICtJ7wzprzNxNOpZ9Vseyz2rxrl/vqHp5uKl9ZKCeH9xObRoxTg5A3VNr/6nn5eWlxMRELV261H7MarVq6dKl6tatm4GVAbVL99gwrX2qjz68t4uevukq3ZbYRPFNAlXP010lZVZtPnxGo2ZtVEZukdGlAoDDOfUdofz8fB04cMD+/NChQ0pJSVFISIiio6M1ceJEjRgxQp07d1ZycrLeeOMNFRQUaNSoUQZWDdQ+Qb5e6h4bpu6xYfZjVqtNh04X6P65m/TjqQLdN2eT5o/tKl8vp/5rAwAqxKmnzy9fvly9e/e+4PiIESM0e/ZsSdKUKVP0yiuvKD09XQkJCXrrrbfUpUuXKn2u2WyW2WyWxWLRvn37mD4Pl3b4dIEGv/ODsgpKdF1cuKYNTZQbO9kDcGIVmT7v1EHIaBX5QQJ12aafsnT39PUqsVg1tmdzPXXjVUaXBACX5BLrCAGoOZ1jQvTKbR0kSe+uPKh5G44YXBEAOAZBCEC53JIQqUf6nduu5ulFO7V6f6bBFQFA1RGEAJTbn/q21OCOkbJYbXrgo83an8HGrQBqN4IQgHIzmUx68Y/tlRQTrLyiMo39cLMK2J4DQC1GELoIs9msuLg4JSUlGV0K4HS8Pdw1bWiiGgWc26vsqYU7xJwLALUVs8Yug1ljwKVt+ilLd7y3TharTc8NaqehXZsaXRIASGLWGIAa0DkmRE9c31qS9I//pGrnsRyDKwKAiiMIAai0MT2aq39cuEosVj3w0WblnC01uiQAqBCCEIBKM5lMevXWeDUJrqe0rLP6y4JtjBcCUKsQhABUSaCvp965p5O83N20JDVD768+ZHRJAFBuBCEAVdahSZCe+cO5bTde/HqPVu47ZXBFAFA+BKGLYPo8UHFDuzbVoIQIlVltGvfhZm09csbokgDgipg+fxlMnwcqpqTMqnvnbNSq/ZkK9vXUgnHd1LKhv9FlAXAxTJ8HYAgvDzdNG5qo+KggnSks1bD3N+h49lmjywKASyIIAXAoP28PzRqZpBYN/HQip0jD3l+vMwUlRpcFABdFEALgcCF+Xpp7bxc1DvTRj6cKNHL2RvYkA+CUCEIAqkVkUD3NHZ2sIF9PbUvL1t++YE8yAM6HIASg2sSG++u9YZ3l7mbSopTjWrDpqNElAcB5CEIAqlVysxBN7N9KkvTslzu1LyPP4IoA4FcEIQDV7oFrW6hHbJiKSq0a/9EWnS2xGF0SAEgiCF0UCyoCjuXmZtLrdySogb+39p/M16QvdxpdEgBIYkHFy2JBRcCxfvgxU0NnrJfVJr1+R7wGd2xidEkA6iAWVATglK5uEaY/9Y2VJP3ti53aeSzH4IoAuDqCEIAa9XCfWF3dIlSFJRbdM2M9YQiAoQhCAGqUu5tJ04YlqmN0kHLOluru6eu0/Wi20WUBcFEEIQA1LsDHU3NHJyuxabByi8p0z4z1SknLNrosAC6IIATAEP4+npozOllJMcHKKyrTsBnrtXp/JqtPA6hRzBq7DGaNAdWvoLhMo2Zv1IZDWZKkhv7e6tW6gXq3bqhrYsMU4ONpcIUAapuK/P4mCF0GQQioGYUlZXp60U59vSNdZ0t/XWzR28NNj13XWqO7N5O7m8nACgHUJgShKjKbzTKbzbJYLNq3bx9BCKghxWUWbTiUpWV7TmnZ3pM6lFkgSUqOCdErt3VQ01A/gysEUBsQhByEO0KAcWw2m+ZvTNM//5uqghKLfL3c9dcbr9I9XaJlMnF3CMClsaAigFrPZDLpzuRoLf5zT3VpFqLCEoueXrRTw2du0Imcs0aXB6COIAgBcGpRIb6aN6arnv1DnLw93LRqf6aue32lFm45ygwzAFVGEALg9NzcTBrdvZn+N6GHEqKClFdUpomfbtPYDzYr52yp0eUBqMUIQgBqjRYN6uuzcd30lwGt5elu0pLUDI37YLNKyqxGlwagliIIAahVPNzdNL53S3027mr5eblr7cHTenrRDrrJAFQKQQhArRQfFaQpd3eSm0n6dNNRTV3xo9ElAaiFCEIAaq3ebRpq0sC2kqSXF+/V/3acMLgiALUNQQhArTbi6hiNvDpGkjThk6165Zs9KiwpM7YoALUGQQhArffMH+L0hw6NVWqxybzsR/WfvFKLd55g3BCAKyIIAaj13N1Mevuujpo2NFGRQfV0LPusxn24Rf/8726jSwPg5AhCF2E2mxUXF6ekpCSjSwFQTiaTSde3a6TvJl6rh3q3lCTNXHNIK/edMrgyAM6MvcYug73GgNrr71/u0uwfflLjQB9980hPBfh4Gl0SgBrCXmMAXN7j17dW01Bfncgp0nP/TTW6HABOiiAEoE7y9fLQq7fFy/TzOkPf78kwuiQAToggBKDOSooJ0b3XNJMkPf7ZDqWkZRtbEACnQxACUKc9NqC1WoXXV2Z+sf449QdN+X6/LFaGRgI4hyAEoE7z8XTXgrFX6w8dGstitenVJft053trdSqv2OjSADgBghCAOi/Q11Nv39VRk2+PV31vD2386YyGvb9e2YUlRpcGwGAEIQAuwWQyaUinJvryoWvU0N9be9LzNGLWRuUXsx0H4MoqtY7Ql19+WeEP6t+/v+rVq1fh9xmJdYSAumlfRp7ueHetzhSWqmvzEM0elSwfT3ejywLgIBX5/V2pIOTmVrEbSSaTSfv371fz5s0r+lGGIggBddeOozm6e/o65RWX6bbEJnrltnijSwLgIDWyoGJ6erqsVmu5vnx9fSv7MQBQLdo3CdS7wxNlMkkLNrPOEOCqKhWERowYUaFurqFDh3JHBYDTubpFmH2doSc/36GcwlKDKwJQ09hr7DLoGgPqvqJSi258a5UOnirQkI6RmnxHgtElAaiiau0aO3PmjLKysiRJp06d0sKFC7Vr167KVQoABvPxdNert8XLzSQt3HpM7yw/oKJSi9FlAaghFQpCM2bMUGJiojp37qypU6dq8ODBWrp0qe68807NmDGjumoEgGrVKTpYY69tIUl6efFe9XpluT5Y+5NKLVaDKwNQ3SrUNdahQwetX79eZ8+eVXR0tA4dOqQGDRooJydH1157rVJSUqqx1JpjNptlNptlsVi0b98+usYAF2C12vTJxjRN+X6/jucUSZK6NAvRO/d0Umh9b4OrA1AR1dY15uHhoXr16ikkJEQtW7ZUgwYNJEmBgYEymUyVr9jJjB8/Xqmpqdq4caPRpQCoIW5uJt3dJVrL/tJL/7ilrep7e2j9oSzdYl6j3SdyjS4PQDWpUBByd3dXUdG5fymtWLHCfjw/P9+xVQGAQbw93DW8W4y+ePBqNQ311dEzZ/XHqT/o3ynHjC4NQDWoUBD67rvv5O197hZxYGCg/XhhYaHee+89x1YGAAaKDffXv8dfo+4tw1RYYtGET1L0ty92MJAaqGOYPn8ZTJ8HYLHa9OZ3+/T2sgOy2aQOTQI1Z1Sygv28jC4NwCXUyMrSv1VaWqq0tDTt3bvXPrUeAOoCdzeTJl7X+lz48fXU9qM5umv6Op3OLza6NAAOUOkglJeXp6lTp+raa69VQECAYmJidNVVV6lBgwZq2rSpxowZw2BjAHVGz1YNtGBcNzX4eef6u6ev16k8whBQ21UqCE2ePFkxMTGaNWuW+vXrp0WLFiklJUX79u3T2rVrNWnSJJWVlem6667T9ddfr/379zu6bgCocS0b+uuT+7uqob+39mbk6dZpP+jgKSaLALVZpcYI3XXXXXr66afVtm3by55XXFysWbNmycvLS6NHj650kUZhjBCAizmUWaDhM9crLeusgnw9NWN4Z3WOCTG6LAA/q8jvbwZLXwZBCMClnMor1n1zNmrb0Rx5uJk06poY/alvrPx9PI0uDXB5NT5YGgBcTQN/b827v6tuat9YZVabpq86pN6vrtCaA5lGlwagAqochF544QXNnDnzguMzZ87USy+9VNXmAcBp+Xp5yHxPJ80alaTmYX7KzC/WuA8268BJxg0BtUWVg9C7776rNm3aXHC8bdu2mjZtWlWbBwCn17t1Q3395x5KjglRXnGZ7p+7STlnS40uC0A5VDkIpaenq3Hjxhccb9CggU6cOFHV5gGgVvD2cNc7QzspItBHBzMLNO6DzTqZW2R0WQCuoMpBKCoqSmvWrLng+Jo1axQREVHV5gGg1gir7633hneWj6eb1h48rb6TV+ij9YfFnBTAeXlUtYExY8boz3/+s0pLS9WnTx9J0tKlS/X444/r0UcfrXKBAFCbtIsM1MIHrtGTC7dr+9Ec/e2LnSqz2DTi6hijSwNwEVWePm+z2fTkk0/qrbfeUklJiWw2m+rVq6cnnnhCzz77rKPqNATT5wFU1i97lL31/QH5eLrpqz/1UIsG9Y0uC3AJhqwjlJ+fr9TUVNWrV0+tWrWy71JfmxGEAFSF1WrT8JkbtPpApuKjgvT5uG7ycGfVEqC61fg6Qu+//766du2qHj16qHPnzkpMTNSMGTMc0TQA1Fpubia9clsHBfh4aFtath5bsE1nSyxGlwXgN6ochJ599llNmDBBAwcO1IIFC7RgwQINHDhQjzzySK3vGgOAqmocWE8v/rGD3EzSopTjGjL1Bx0+XWB0WQB+VuWusQYNGuitt97SXXfddd7xefPm6eGHH1ZmZu1dZZWuMQCOsvbH03p43hZl5peoob+3Ph3bTTFhfkaXBdRJNdo1Vlpaqs6dO19wPDExUWVlZVVtHgDqhG4tQvXfh3uodbi/TuYV654Z67UnPVdWK1PrASNVOQgNGzZMU6dOveD4e++9p3vuuaeqzRvCbDYrLi5OSUlJRpcCoA5pFOijD+/rouYN/HQs+6yuf2OV2v39Gz3+2TYVlTJ2CDBClbvGHn74Yc2dO1dRUVHq2rWrJGn9+vU6cuSIhg8fLk/PX3dinjx5ctWqrWF0jQGoDuk5Rfrz/K3afPiMSi3n/gru1jxU7w1PZPd6wAFqdPp87969y3WeyWTS999/X5WPqnEEIQDVqdRi1er9mXp43lblF5epU3SQPh3LFHugqgxZR6guIggBqAk7j+XorunrlFdUpkkD4zTqmmZGlwTUajW+jhAAoPLaRQbqyRvaSJJeW7JPGWzWCtSYSu81Nnr06HKdN3PmzMp+BAC4jLuSovXZ5qPaeiRbw9/foP5x4RrSKVLN2ZYDqFaV7hpzc3NT06ZN1bFjx8vurPzFF19Uujij0TUGoCbtOp6jW6eu1dmfZ5DV9/bQ9OGd1a1FqMGVAbVLjYwRGj9+vObNm6emTZtq1KhRGjp0qEJCQipVsLMiCAGoaUfPFGrlvkwt2JymrUey5eXhpjfvSNAN7RsbXRpQa9TIGCGz2awTJ07o8ccf13/+8x9FRUXp9ttv1zfffHPZO0QAgEtrEuyru7tEa96YruofF66SMqse+GiL3lq6n79bgWrgsFljhw8f1uzZszV37lyVlZVp165dql+/dvdtc0cIgJHKLFY999Vuzf7hJ0lSu8gATezfSn3ahBtbGODkDJk15ubmJpPJJJvNJouFFVIBoKo83N3095vb6sUh7eXr5a6dx3I1evYmfboxzejSgDqjSkGouLhY8+bNU//+/dWqVSvt2LFDU6ZM0ZEjR2r93SAAcBZ3Jkdr9RN9dFdytCTp2S936sDJPIOrAuqGSk+ff/DBB/XJJ58oKipKo0eP1rx58xQWFubI2gAAPwvx89Lzg9opLatQqw9kauiMDbqhfSPd16O5IoPqGV0eUGtVafp8dHS0OnbsKJPJdMnzFi5cWOnijMYYIQDO5mRekQZNWaPjOecWXWwW5qdv/txTXh6sjwv8oiK/vyt9R2j48OGXDUAAAMdr6O+jbx7pqVX7MzXpy106lFmgOT/8pDE9mxtdGlArsdfYZXBHCIAz+3RTmh7/bLv8vT303z91V9NQP6NLApwCe40BgAu4tVMTtY8MVF5xmfq/vlKT/r1TP/yYKauVf98C5eWwILR+/XpHNQUAKAc3N5PeuaeTrm4RqpIyq+asPay7p6/XyNkbVVBcZnR5QK3gsK6x6OhoHTlyxBFNOQ26xgDUBjabTSv2ndJ/tp3QVzuOq6jUqg5NAvXRfV3k7+NpdHlAjau2wdK33377RY/bbDZlZWVVpCkAgIOYTCb1at1QvVo31LBuTTV69kZtP5qjRz/dpmlDE+XmxsQW4FIqdEcoJCREH3zwwQWLJdpsNt1xxx3KyMhweIFG4o4QgNpoW1q2bpu2ViUWq+5KjtKj17VWWH1vo8sCaky13RHq1auX/P391bNnzwte69ChQ8WqBABUi/ioIP1zUFs98fkOzduQpi9Tjuu5we00uGMTo0sDnA7T5y+DO0IAarPle0/qtSX7tONYjiTpjs5R+vvNbVXPy93gyoDqVWPT59PT06vydgBANerVuqEWjb9Gj/RrJZNJmr8pTYPMa3TgZL7RpQFOo0pB6LrrrnNUHQCAauDuZtKEfrH66N4uCqvvrb0Zebplymp9ue240aUBTqFKQYheNQCoHa5uGab/Teiubs1DVVBi0Z/mbdWNb67S9W+sVOrxXKPLAwxTpSDEXmMAUHs09PfRh/d10cN9WspkklJP5GpPep4eXbBNpRar0eUBhmCLDQBwIe5uJj16XWt99XAPvXlngoJ8PbX7RK6mLf/R6NIAQ1R693kAQO0VFxGguIgAlZRZ9ZfPtuu1b/epnpe77uvBLvZwLVW6I+TuzhRMAKjNbk1sonHXtpAkPffVbv31ix0qLrMYXBVQc6oUhLZu3eqoOgAABjCZTHri+tZ64vo2Mpmkj9cf0V3vrdPJ3CKjSwNqRJ0eI5Sdna3OnTsrISFB7dq10/Tp040uCQCcjslk0gO9WmjmyCQF+Hhoy5Fs/eHt1Vq+96QsVmYHo3qUWqzacChLJWXGDtSv0ytLWywWFRcXy9fXVwUFBWrXrp02bdqk0NDQcr2flaUBuJqfMgt0/webtC/j3KKLkUH19P7IzmrTiL8D4Vh/+2KHPlp/RHd0jtJLtzp2m64aW1lakjZu3Ki+ffuqQ4cOGjJkiP7xj3/oyy+/1JEjR6radJW5u7vL19dXklRcXCybzcbaRwBwGTFhfvriwWs0tGu0Anw8dCz7rO58b53WHMg0ujTUMR+tP5cT5m9KM7SOKgehYcOGyd3dXffff7+aNWumFStWaNSoUYqJiSn3nZdLWblypQYOHKiIiAiZTCYtWrTognPMZrNiYmLk4+OjLl26aMOGDee9np2drfj4eDVp0kR/+ctfFBYWVqWaAKCu8/P20HOD2mvVE30UHxWk7MJS3TNjvcZ9sFn7MvKMLg9wqCpPn09LS9NXX32lFi1anHf88OHDSklJqVLbBQUFio+P1+jRozVkyJALXp8/f74mTpyoadOmqUuXLnrjjTc0YMAA7d27Vw0bNpQkBQUFadu2bcrIyNCQIUN06623Kjw8vEp1AYArCKznqXljuujFr/do7trDWrwrXesPndZXf+qhiKB6RpcHOESV7whdc801Onr06AXHmzZtqltuuaVKbd9www167rnnNHjw4Iu+PnnyZI0ZM0ajRo1SXFycpk2bJl9fX82cOfOCc8PDwxUfH69Vq1Zd8vOKi4uVm5t73hcAuDJfLw/945Z2WvJIT8U1DtCZwlL9ceoPenrRDqXnMLMMtV+lgtCQIUP097//XV988YXGjRunf/7znzpz5oyja7uskpISbd68Wf369bMfc3NzU79+/bR27VpJUkZGhvLyzt3GzcnJ0cqVK9W6detLtvnCCy8oMDDQ/hUVFVW93wQA1BKtwv01bWiiQv28dCKnSB+uO6Lery6XedkBlbE9B2qxSnWNtWjRQmvWrNE777yjzMxzA+hatWqlW265RV27dlXHjh3Vvn17eXl5ObTY38rMzJTFYrmgmys8PFx79uyRdK577v7777cPkn744YfVvn37S7b51FNPaeLEifbnubm5hCEA+Fl0qK++f7SX1h86rXdXHtTmw2f0yjd7lZ5TpH8Oamd0eUClVCoIvfLKK/bHx44dU0pKiv3rpZde0sGDB+Xh4aHWrVtr+/btDiu2opKTkys0Tsnb21ve3t7VVxAA1HKBvp66rm0j9Y8L16eb0vTE5zv0wbrD8vVy1yP9W8nHkx0HULtUebB0ZGSkIiMjddNNN9mP5efnKyUlRdu2batq85cUFhYmd3d3ZWRknHc8IyNDjRo1qrbPBQCcW4TxjqRoZeQWa/K3+/TuyoNad/C0zPd0UpNgX6PLQy3i4WYy9PMrNUboSmsE1a9fX927d9f48eMlnbtr5GheXl5KTEzU0qVL7cesVquWLl2qbt26OfzzAAAX+lPfWL03LFHBvp7adjRHPV9epslL9hpdFmoR99oYhJKSkjR27Fht3Ljxkufk5ORo+vTpateunT7//PNKFffLnaVfurcOHTqklJQUexCbOHGipk+frjlz5mj37t164IEHVFBQoFGjRlXq8wAAFXdd20Za+OA16to8RFab9Nb3BzTl+/0qZRA1LuG3ixsbfUeoUl1jqampev7559W/f3/5+PgoMTFRERER8vHx0ZkzZ5Samqpdu3apU6dOevnll3XjjTdWqrhNmzapd+/e9ue/DGQeMWKEZs+erTvuuEOnTp3Ss88+q/T0dCUkJGjx4sVVXifIbDbLbDbLYmEHZgAoj2Zhfvrk/m56a+l+Tf52n15dsk9LUjP095vbqn1koDzd6/TWlqigjT/9OtPc6P0eqrTX2NmzZ/XVV19p9erVOnz4sM6ePauwsDB17NhRAwYMULt2tXsWAXuNAUDF2Gw2fbzhiF76eo9yi8okSVEh9TRjeJJaN/I3uDo4i1X7T2nY++d2gjCZpEMv3HSFd1RMRX5/1+lNV6uKIAQAlZN6PFePzE/R3p+35PDzctebd3ZUvzhW9oe0fO9JjZz16/Can140LghxrxIA4HBxEQH65pGe2vpMf3VtHqKCEovGfLBJT36+nf3KoD3pznMNVCoIbd++XVZr+QfB7dq1S2VlZZX5KABALRbs56UP7u2ie7pEy2aTPtmYphvfXKVXvtmjolLGYbqqF7/eY3QJdpUKQh07dtTp06fLfX63bt2uOOUeAFA3ebq76fnB7fXhvV3U76qGKrPaZF72o8Z+sJntOaA2Bo8dq9SsMZvNpmeeeUa+vuVbNKukpKQyH2MYZo0BgON1jw1T99gwLd55Qo/M36YV+05p4JQ1mtA3VgPahstkMnYaNYzRtXmooZ9fqcHSvXr1qvAF+/HHH6tx48YV/ShDMVgaAKrHd6kZmvDJVhWUnPsH5zUtQ/XcoPZqFuZncGWoCTFPfmV/PKJbU/3fLY6dZV6R39+VuiO0fPnyyrwNAABJUr+4cK1+oo9mrD6oGasOac2B07rxzVV67fZ4XRcXLg/WHUIN4UoDABgi2M9LfxnQRkse6ankmBCdLbXowY+26OYpa3Qi56zR5cFFEIQAAIZqGuqnufcm686kKHl7uCn1RK5ufHOVXl68R1kFtWuMKa4s52zpec+NXsyQIAQAMJyPp7te/GMHff9YL7UKr68zhaV6Z/mPunXaDzp6ptDo8uBALy92nqnzEkEIAOBEIoPq6bMHrtbTN12lBv7eOniqQEPe+UHfpmbIajX63gEc4bPNR40u4TwEoYswm82Ki4tTUlKS0aUAgMsJ8PHUfT2a68uHrlHrcH+dzCvWmLmbdPeMdcrMLza6PFRRp+hgo0s4j0OCUGlpqdLS0rR3715lZWU5oklDjR8/Xqmpqdq4ceOVTwYAVIvGgfW08MGrdW/3Zqrn6a51B7PU+5Xl+r//7GIhxlosIqjeec+N3vG00kEoLy9PU6dO1bXXXquAgADFxMToqquuUoMGDdS0aVONGTOGIAEAqBI/bw8984c4ffnQNWrewE95xWWateYn3f/BZv13+3Gl5xQZXSIqKL/43GDpUD8vgys5p1JBaPLkyYqJidGsWbPUr18/LVq0SCkpKdq3b5/Wrl2rSZMmqaysTNddd52uv/567d+/39F1AwBcSGy4vxZP6KnnBrWTh5tJ3+85qYc+3qrery7XzmM5RpeHCigoPreIpp93pZYydLhKVbFx40atXLlSbdu2vejrycnJGj16tKZNm6ZZs2Zp1apVio2NrVKhAADX5uXhpqFdm6ptRIDeX31Iu47n6lBmgcZ+sFkzRyaptcF7VqF8Vh/IlCTV/zkI2QyeQF+pIDRv3jz745MnT6phw4YXPc/b21vjxo2rXGUAAFxEx+hgTbk7WDmFpbrZvFqHTxfqhjdXavQ1zTSmZ3OFB/gYXSIuYVtatv1xfR/nuCNU5cHSt9566yU3Jy0rK6tq8wAAXFSgr6fmjEpWz1YNZLVJM1YfUo+Xl+n7PRlGl4ZL2P6bbsygep4GVvKrKgehoKAg/elPf7rg+OnTp9WvX7+qNg8AwCXFhPlp7uhkzRqZpI7RQSops2r07E16auEO5RWVXrkB1KjfLqbYNNTXwEp+VeUgNHfuXH377beaOXOm/dju3buVnJwsP7/auYsw6wgBQO3Su01Dzb+/m/pddW6oxrwNR9T+70s0/uMtKiq9eK8Fal5e0a89Rb5eP48Rqq3T538RFBSkzz//XH/5y1+0YcMGffPNN+rWrZsGDRqk//znP46oscaxjhAA1D5eHm6aMSJJn9zfVcG+57pdvtp+QsPeX880eyewNz3vvOcmk0GF/E6lRioNGTJECQkJ9q/27dtrypQpuvHGG1VUVKS3335bo0aNcnStAABcUdfmofrmzz315bbjenPpfm386Yy6vbhUY3u20GPXtZKHO5sqGOH0b1YFn9DXeWaSV+pqaNGihVatWqX77rtPMTExCg0N1fTp02Wz2XT33XerU6dOKi2lbxYAYIyGAT66r0dzfTq2m2Ib1pfNJk1b8aNa/u1rPfjRZlamNtgN7RvZHxu9g1yl7gi98sor9sfHjh1TSkqKUlJSFBoaqmXLlun999+Xh4eH2rRpo23btjmsWAAAKuKqxgH6duK1+t+OE3ris+3KKy7T/3akS9qq1+9IkLeHu9Eluowzhb/eIGnTKEDf7HSO2X1VnsQfGRmpyMhI3XTTTfZj+fn5SklJIQQBAJzCje0bq31koG5/d61O5BTpfzvS9cOPS/XKrfHqHxdudHkuYfzHWyRJN/7mbpAzqFTX2JEjRy77ev369dW9e3eNHz9e0rm7RgAAGCkqxFcr/tJbE/rGytPdpOzCUo2Zu0nP/TdVNqOnLtVxv/359mnjXMGzUkEoKSlJY8eOveysqpycHE2fPl3t2rXT559/XukCAQBwFC8PNz3Sv5XWPNFHN8dHSDq3EOPo2RvPG8wLx/ptt9jA+MbnvWZ0Bq1U11hqaqqef/559e/fXz4+PkpMTFRERIR8fHx05swZpaamateuXerUqZNefvll3XjjjY6uGwCASmsY4KO37uqopJhg/fOr3Vq295R6vrxM9/Vorof6tJQnM8sc6pflC0L8vOzjspxl+nyl/qRDQ0M1efJknThxQmazWbGxscrMzLTvMn/PPfdo8+bNWrt2LSEIAOC0hnWL0cIHrlbbiAAVlFj05tL9uv6Nlfp3yjFZrXSXOcqOY9mSpKjgesYWchFVGix98uRJeXl56e6771ZycrKjajKc2WyW2Wy+5B5qAIC6o11koP77cHd9ue24nvh8u348VaAJn6TofztO6MUhHRTs52V0ibXet6knJekSA9ONDZyVvvc3b948tWrVSrfccou6deumzp0769SpU46szTCsLA0ArsVkMumWhEh9Nu5q3ZrYRB5uJn2zK0NJz3+nD9cdZjB1FeQUluq73eemyndrEWo/7iQ9Y5UPQv/3f/+nu+++W3v27NGSJUskSU8++aTDCgMAoKa1iwzUq7fF673hiQqr760yq01PL9qpkbM26mQe23RUxta0M/bHbSMCDazk4iodhA4ePKhJkyapVatW6tu3rz788EN98sknjqwNAABD9GkTro1/66tH+rWSJK3Yd0qDzT9o0dZj3B2qoF3HcyVJSTHB8vF0vgUsKx2EysrK5Ovra3/epk0bWa1WpaenO6QwAACMZDKZNKFfrGYM76yw+t46ln1Wf56folGzN+pEzlmjy6s1dh7LkST1btPwoq8bnSurND9wzpw5+uGHH5Sfny9J8vDwUGFhoUMKAwDAGfSLC9fKx3vZNwpdvveU+ry6Qi8t3sPdoSuw2Wz6eue5GyRRwb7nvVarp89LUo8ePfTcc8+pe/fuCgoKUmxsrIqKivT+++9r2bJlysvLc2SdAAAYxtfLQ4/0b6V5Y7qqTSN/nS21aOryH/XUwh3cHbqM1BO59sfJzUIMrOTSKh2EVqxYoZycHO3du1cffvihBg8erGuvvVZTp05V3759FRwcrKuuusqRtQIAYKhuLUL19YQeevKGNpKkTzamqe9rK/T1jhMGV+ac/vGfVElSz1YNFB7gY3A1F1flTVdjY2MVGxurO++8037s0KFD2rRpk7Zu3VrV5gEAcComk0njrm2hVuH19a//7dGBk/l64KMt6t4yTC8Maa+oEN8rN+ICcgpLtf5QliTp9s5NLnme0b2L1bKGeLNmzXTbbbfpX//6V3U0DwCA4fq0CdfiCT00pkczuZmk1QcyNfidH1iV+mefbkqzP+7aPPSC101OMkiIzVQAAKgkD3c3/e2mOC17rJdaNqyvzPxiTfgkRXdOX+fy6w7tyzg3VrhpqK/C6nsbXM2lEYQAAKiipqF++nRsN93XvZm83N204VCWkp9fqqnLfzS6NEPYbDZtOnxuIcXHB7S5/Lm1dYuNusxsNisuLk5JSUlGlwIAqCVC/Lz09B/i9PWfeygy6Nzmoi8t3qOxH2zS8WzXmln24EdbdCizQD6eburVuoHR5VwWQegi2GsMAFBZLRrU18IHr1Zsw/qSpG92ZWjA6yu1an/d2I/zSopKLfa1gzpEBsnPu8rzsqoVQQgAAAcLD/DRtxOv1Uf3dVHrcH/lFZdp2PsbNGbuJh0+XWB0edXq7e/32x+/elu8gZWUD0EIAIBqck3LMH3+4NW6u0u03EzSt6kZuumt1fpkw5E6ObOsqNSiGasOSZKGd2uq6NArLyVQJ6fPAwCAc+p7e+hfg9trySM9lRQTrPziMj25cIcGT/3Bvg9XXfHcV6kqLrNKkn3D2ktxktnzBCEAAGpCy4b++nhMVz1901Wq7+2hbWnZ+sPbq3Xne2tr7VT7X/Zas1ptuv6Nlfpw3RFJ0lM3tFGwn5eRpZWbc49gAgCgDvF0d9N9PZrr5vgIPffVbn257bjWHTw31f6ft7TVrYlRquflbnSZV7Ri3ym9vHiPdp/IVZ824fpud4b9tZ6tGuj+ns3L3ZbRHYTcEQIAoIY1DPDRW3d11PsjOtuPPfPvXRrwxkp9m5ohixOPH/rnf1M1YuYG7TqeK6tN54UgSZoxvHO5Vo02yTn6xghCAAAYpO9V4Vr7VB8NSohQkK+njmQVaszcTbrxzVVavPOEyixWo0s8z46jOXp/9aGLvjbqmhj99OJN8vKoXdGCrjEAAAzUOLCe3rizo/KKSjX52336bNNR7c3I07gPtyjY11Ojr2mm25OiyrV7u81mU35xmep7ezh8L69le0/qkfkp9ufb/36dfDzclZFbpCbB9Zxm77CKIggBAOAE/H08NWlgW03oG6t3Vx7U+6sP6UxhqV77dp/e/v6AerZqoKtbhKp/XPgFweN49lkt33tK8zelaVtatvx9PPS3G6/SncnRDqltya503f/BZklShyaBmjGiswJ8PCVJUSFXniJ/OUZPnzfZbEaX4Lxyc3MVGBionJwcBQQEGF0OAMCFlFqsmr8xTQs2H9W2tOwLXg/29ZTJZFKpxaq8orKLtpHYNFh/vbGNEpuGVLqOpbszdO+cTfb2Prqvi3w8qz6ge9qKH/Xi13v0x05N9Nrtjl14sSK/v7kjBACAE/J0d9PQrk01tGtTbfwpS19sPaaFW46qqPTcuKEzhaX2c91MUqfoYPVu01C9WzfU8JnrlZlfos2Hz+iPU9fq+raN9MQNbdQszK/cn2+12nT3jHVadzDLfmzq0E4OCUHOhCAEAICTS4oJUVJMiJ4f1E4ZucVavvekgv285OvlrgAfT8WE+inQ19N+/sa/9dOGQ1n6cP0R/WfbcS3ela7Fu9J1VeMA3dS+kQZ3amLfGPZithw5o5cX7zkvBG34a1819L/yOKWKMnr3eYIQAAC1hMlkUqNAnyuO/TGZTOrSPFRdmodq5NUxmvL9fi3be0q7T+Rq94lcvbpkn4J8PdUxKkhdm4fK3c2ktKxC5Rdb9OW2Yyq1/BpOerZqoNkjk+Tm5tjB0M4ytJogBABAHZbYNFizRiVr94lcrd6fqaV7MrT+UJayC0u1bO8pLdt76qLvS4oJ1r8Gt1dsuH8NV1yzCEIXYTabZTabZbFYjC4FAACHuKpxgK5qHKAxPZvrZF6RPvt5EHZxmVXZhaUKD/BWTKifokN91TrcX4lNg2vtlPiKIAhdxPjx4zV+/Hj7qHMAAOqShv4+erBXS6PLOIfd5wEAgKtxlptNBCEAAOCyCEIAAMAwRq/qTBACAAA1jt3nAQAADEYQAgAALosgBAAADGP03u8EIQAAUOOYPg8AAGAwghAAADAM0+cBAAAMQhACAAAuiyAEAABcFkEIAAAYxuDZ8wQhAABQ80xOMn+eIAQAAFwWQQgAABiG6fMAAMDlOEfHGEEIAAC4MIIQAABwWQShizCbzYqLi1NSUpLRpQAAUKex+7wTGj9+vFJTU7Vx40ajSwEAoE5yktnzBCEAAOC6CEIAAMAwTJ8HAAAux0l6xghCAADAdRGEAACAyyIIAQAA47D7PAAAcDXsPg8AAGAwghAAADCMzeC+MYIQAACocU7SM0YQAgAArosgBAAAXBZBCAAAGMbgzecJQgAAoOY5yRAhghAAAHBdBCEAAGAYusYAAIDrcZL58wQhAADgsghCAADAZRGEAACAYdhiAwAAuBznGCFEEAIAAC6MIAQAAAzD9HkAAOBynGT2PEEIAAC4LoIQAABwWQQhAABgGIOHCBGEAABAzTM5yQR6ghAAAHBZdToIpaWlqVevXoqLi1OHDh20YMECo0sCAABOxMPoAqqTh4eH3njjDSUkJCg9PV2JiYm68cYb5efnZ3RpAABAxq8jVKeDUOPGjdW4cWNJUqNGjRQWFqasrCyCEAAABmMdoXJYuXKlBg4cqIiICJlMJi1atOiCc8xms2JiYuTj46MuXbpow4YNF21r8+bNslgsioqKquaqAQBAbeHUQaigoEDx8fEym80XfX3+/PmaOHGiJk2apC1btig+Pl4DBgzQyZMnzzsvKytLw4cP13vvvVcTZQMAgHIztm/MqbvGbrjhBt1www2XfH3y5MkaM2aMRo0aJUmaNm2avvrqK82cOVNPPvmkJKm4uFiDBg3Sk08+qauvvvqyn1dcXKzi4mL789zcXAd8FwAA4PecpGfMue8IXU5JSYk2b96sfv362Y+5ubmpX79+Wrt2rSTJZrNp5MiR6tOnj4YNG3bFNl944QUFBgbav+hGAwCgbqu1QSgzM1MWi0Xh4eHnHQ8PD1d6erokac2aNZo/f74WLVqkhIQEJSQkaMeOHZds86mnnlJOTo79Ky0trVq/BwAAYCyn7hqrqu7du8tqtZb7fG9vb3l7e1djRQAA4LeMnj5fa+8IhYWFyd3dXRkZGecdz8jIUKNGjQyqCgAAlAfT56vIy8tLiYmJWrp0qf2Y1WrV0qVL1a1bNwMrAwAAtYVTd43l5+frwIED9ueHDh1SSkqKQkJCFB0drYkTJ2rEiBHq3LmzkpOT9cYbb6igoMA+iwwAADg3o3efd+ogtGnTJvXu3dv+fOLEiZKkESNGaPbs2brjjjt06tQpPfvss0pPT1dCQoIWL158wQDqijKbzTKbzbJYLFVqBwAAXJyz7D5vstmMHqbkvHJzcxUYGKicnBwFBAQYXQ4AAHXGpxvT9Pjn29WnTUPNHJnk0LYr8vu71o4RAgAAqCqCEAAAMIzRHVMEIQAAUPOcY4gQQQgAALgugtBFmM1mxcXFKSnJsYO3AADA+YyesUUQuojx48crNTVVGzduNLoUAADqJCfpGSMIAQAA10UQAgAALosgBAAADGP0ss4EIQAAUONMTrL9PEEIAAC4LILQRTB9HgCAmsH0eSfE9HkAAKqXc3SMEYQAAIALIwgBAACXRRACAACGYfd5AADgcpxk9jxBCAAAuC6CEAAAcFkEoYtgHSEAAKoXXWNOjHWEAABwDQQhAADgsghCAADAMOw+DwAAXI7JSTbZIAgBAACXRRACAAA1LsTPS0kxwWoV7m9oHR6GfjoAAHBJPVs1UM9WDYwugztCAADAdRGEAACAyyIIXQQrSwMA4BpMNpvRM/idV25urgIDA5WTk6OAgACjywEAAOVQkd/f3BECAAAuiyAEAABcFkEIAAC4LIIQAABwWQQhAADgsghCAADAZRGEAACAyyIIAQAAl0UQAgAALosgBAAAXJaH0QU4I7PZLLPZrLKyMknnluoGAAC1wy+/t8uzixh7jV3G0aNHFRUVZXQZAACgEtLS0tSkSZPLnkMQugyr1arjx4/L399fJpNJkpSUlKSNGzde8b2OOC83N1dRUVFKS0urc5u+lvfnUxs/3xFtV6WNir7Xkdf0lc7hmq59n++odivbTnVdz+U911X/jpaMvaar+tk2m015eXmKiIiQm9vlRwHRNXYZbm5uFyRJd3f3cl3wjjwvICCgzv1PVt6fT238fEe0XZU2KvpeR16r5W2La7r2fL6j2q1sO9V1PZf3XFf9O1oy9pp2xGcHBgaW6zwGS1fQ+PHjDTmvrjH6+67Oz3dE21Vpo6LvdeS1avSfq5GM/t6r6/Md1W5l26mu67m85xr952okI7/3mvxsusacWG5urgIDA5WTk1Mn/7UB18M1jbqE67lu4I6QE/P29takSZPk7e1tdCmAQ3BNoy7heq4buCMEAABcFneEAACAyyIIAQAAl0UQAgAALosgBAAAXBZBqA5IS0tTr169FBcXpw4dOmjBggVGlwRU2eDBgxUcHKxbb73V6FKASvnvf/+r1q1bKzY2VjNmzDC6HFwCs8bqgBMnTigjI0MJCQlKT09XYmKi9u3bJz8/P6NLAypt+fLlysvL05w5c/TZZ58ZXQ5QIWVlZYqLi9OyZcsUGBioxMRE/fDDDwoNDTW6NPwOd4TqgMaNGyshIUGS1KhRI4WFhSkrK8vYooAq6tWrl/z9/Y0uA6iUDRs2qG3btoqMjFT9+vV1ww03aMmSJUaXhYsgCNWAlStXauDAgYqIiJDJZNKiRYsuOMdsNismJkY+Pj7q0qWLNmzYUKnP2rx5sywWi6KioqpYNXBpNXlNA0ao6jV+/PhxRUZG2p9HRkbq2LFjNVE6KoggVAMKCgoUHx8vs9l80dfnz5+viRMnatKkSdqyZYvi4+M1YMAAnTx50n5OQkKC2rVrd8HX8ePH7edkZWVp+PDheu+996r9e4Jrq6lrGjCKI65x1BI21ChJti+++OK8Y8nJybbx48fbn1ssFltERITthRdeKHe7RUVFth49etjmzp3rqFKBcqmua9pms9mWLVtm++Mf/+iIMoFKq8w1vmbNGtugQYPsr0+YMMH20Ucf1Ui9qBjuCBmspKREmzdvVr9+/ezH3Nzc1K9fP61du7ZcbdhsNo0cOVJ9+vTRsGHDqqtUoFwccU0Dzqw813hycrJ27typY8eOKT8/X19//bUGDBhgVMm4DIKQwTIzM2WxWBQeHn7e8fDwcKWnp5erjTVr1mj+/PlatGiREhISlJCQoB07dlRHucAVOeKalqR+/frptttu0//+9z81adKEEAWnUZ5r3MPDQ6+99pp69+6thIQEPfroo8wYc1IeRheAquvevbusVqvRZQAO9d133xldAlAlN998s26++Wajy8AVcEfIYGFhYXJ3d1dGRsZ5xzMyMtSoUSODqgIqj2sadR3XeN1CEDKYl5eXEhMTtXTpUvsxq9WqpUuXqlu3bgZWBlQO1zTqOq7xuoWusRqQn5+vAwcO2J8fOnRIKSkpCgkJUXR0tCZOnKgRI0aoc+fOSk5O1htvvKGCggKNGjXKwKqBS+OaRl3HNe5CjJ625gqWLVtmk3TB14gRI+znvP3227bo6Gibl5eXLTk52bZu3TrjCgaugGsadR3XuOtgrzEAAOCyGCMEAABcFkEIAAC4LIIQAABwWQQhAADgsghCAADAZRGEAACAyyIIAQAAl0UQAgAALosgBAAAXBZBCAAAuCyCEADUgMGDBys4OFi33nqr0aUA+A2CEADUgAkTJmju3LlGlwHgdwhCAGrUY489pkGDBpXr3F69eslkMslkMiklJaVSbTiLXr16yd/f/6KvjRw50v59Llq0qGYLA1wcQQhAjUpJSVFCQkK5zx8zZoxOnDihdu3anddGhw4dLvmeX4LFuHHjLnht/PjxMplMGjlyZEXKrlZvvvmmTpw4YXQZgEsiCAGoUdu2batQEPL19VWjRo3k4eFxXhvx8fGXfV9UVJQ++eQTnT171n6sqKhIH3/8saKjoytc95UkJCSoXbt2F3wdP378iu8NDAxUo0aNHF4TgCsjCAGoMUePHlVmZqY9CGVnZ2vgwIHq3r270tPTK9SGJPXv31++vr5q3bq11q9ff955nTp1UlRUlBYuXGg/tnDhQkVHR6tjx47nndurVy899NBDeuihhxQYGKiwsDA988wzstls9nOsVqtefvlltWzZUt7e3oqOjtbzzz9vfz0lJUU7d+684CsiIqJCPyMANYsgBKDGpKSkKCgoSDExMdqxY4eSkpIUGRmpZcuWlfuOyC9jhcxms/76179q27Ztio6O1pNPPnnBuaNHj9asWbPsz2fOnKlRo0ZdtN05c+bIw8NDGzZs0JtvvqnJkydrxowZ9tefeuopvfjii3rmmWeUmpqqjz/+WOHh4RX47gE4I48rnwIAjpGSkqL4+Hh9/PHHeuihh/TSSy9pzJgxFW4jJCREn376qcLCwiRJN998s959990Lzh06dKieeuopHT58WJK0Zs0affLJJ1q+fPkF50ZFRen111+XyWRS69attWPHDr3++usaM2aM8vLy9Oabb2rKlCkaMWKEJKlFixbq3r17uevu16+ftm3bpoKCAjVp0kQLFixQt27dKvS9A3A8ghCAGpOSkqLt27froYce0ldffVWpIJCSkqJbbrnFHoIk6dChQ2rZsuUF5zZo0EA33XSTZs+eLZvNpptuuum89/1W165dZTKZ7M+7deum1157TRaLRbt371ZxcbH69u1b4Xp/8d1331X6vQCqD11jAGpMSkqKhgwZoqKiImVnZ1e6ja5du15w7FIDsEePHq3Zs2drzpw5Gj16dKU+s169epV6HwDnRxACUCPy8vJ08OBBjR8/XlOmTNGdd96pXbt2VaqN3w92vlwQuv7661VSUqLS0lINGDDgkm3/frD1unXrFBsbK3d3d8XGxqpevXpaunRpheoF4PzoGgNQI7Zt2yZ3d3fFxcWpY8eO2rlzpwYOHKgNGzZcsrvqUm20b9/efuzw4cM6c+bMJYOQu7u7du/ebX98KUeOHNHEiRM1duxYbdmyRW+//bZee+01SZKPj4+eeOIJPf744/Ly8tI111yjU6dOadeuXbr33nvL+RMA4IwIQgBqREpKitq0aSNvb29J0iuvvKLdu3dryJAh+u677+Tl5VWuNlq3bi0fHx/7sa1bt9pnol1KQEDAFdsePny4zp49q+TkZLm7u2vChAm6//777a8/88wz8vDw0LPPPqvjx4+rcePGF12wEUDtYrL9dqEMAHAivXr1UkJCgt5444068TlXYjKZ9MUXX9S67UOA2owxQgCc2jvvvKP69etrx44dRpdSbcaNG6f69esbXQbgkrgjBMBpHTt2zL5FRnR0dLm6zyrD6DtCJ0+eVG5uriSpcePG8vPzM6QOwBURhAAAgMuiawwAALgsghAAAHBZBCEAAOCyCEIAAMBlEYQAAIDLIggBAACXRRACAAAuiyAEAABcFkEIAAC4LIIQAABwWQQhAADgsghCAADAZRGEAACAy/p/JzMux1Rcs0EAAAAASUVORK5CYII=\n"},"metadata":{}}]},{"metadata":{},"cell_type":"markdown","source":"### Now lets compute the power spectrum of halos with masses above 1e14 in redshift-space"},{"metadata":{"trusted":true},"cell_type":"code","source":"snapdir = '/home/jovyan/Data/Halos/FoF/s8_p/397/' #folder hosting the catalogue\nsnapnum = 4                                       #number of the catalog (4-->z=0, 3-->z=0.5, 2-->z=1, 1-->z=2, 0-->z=3)","execution_count":7,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Lets read the halo catalog"},{"metadata":{"trusted":true},"cell_type":"code","source":"# read the halo catalogue\nFoF = readfof.FoF_catalog(snapdir, snapnum, long_ids=False,\n                          swap=False, SFR=False, read_IDs=False)\n\n# get the properties of the halos\npos_h  = FoF.GroupPos/1e3            #Halo positions in Mpc/h\nvel_h  = FoF.GroupVel*(1.0+redshift) #Halo peculiar velocities in km/s\nmass_h = FoF.GroupMass*1e10          #Halo masses in Msun/h\nNp_h   = FoF.GroupLen                #Number of CDM particles in the halo. Even in simulations with massive neutrinos, this will be just the number of CDM particles","execution_count":8,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Lets move halos to redshift-space along the z axis:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# move halos to redshift-space. After this call, pos_h will contain the\n# positions of the halos in redshift-space\naxis = 2 #axis along which to displace halos\nRSL.pos_redshift_space(pos_h, vel_h, BoxSize, Hubble, redshift, axis)","execution_count":9,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Lets now select all halos with masses above 1e14 Msun/h"},{"metadata":{"trusted":true},"cell_type":"code","source":"indexes = np.where(mass_h>1e14)[0]\npos_h  = pos_h[indexes]\nvel_h  = vel_h[indexes]\nmass_h = mass_h[indexes]\nNp_h   = Np_h[indexes]\n\nprint('%.3e < Mass M < %.3e Msun/h'%(np.min(mass_h), np.max(mass_h)))\nprint('%d halos with masses above 1e14 Msun/h'%pos_h.shape[0])","execution_count":10,"outputs":[{"output_type":"stream","text":"1.005e+14 < Mass M < 4.589e+15 Msun/h\n39096 halos with masses above 1e14 Msun/h\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"Now lets construct the density field of these halos"},{"metadata":{"trusted":true},"cell_type":"code","source":"# define the matrix hosting the density field\ndelta_h = np.zeros((grid,grid,grid), dtype=np.float32)\n\n# construct 3D density field\nMASL.MA(pos_h, delta_h, BoxSize, MAS, verbose=verbose)\n\n# compute the overdensity field\ndelta_h /= np.mean(delta_h, dtype=np.float64);  delta_h -= 1.0\n\n# print some information\nprint('%.3f < delta < %.3f'%(np.min(delta_h), np.max(delta_h)))\nprint('< delta > = %.3f'%np.mean(delta_h))","execution_count":11,"outputs":[{"output_type":"stream","text":"\nUsing CIC mass assignment scheme\nTime taken = 0.248 seconds\n\n-1.000 < delta < 5290.054\n< delta > = 0.000\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"We can compute the power spectrum now"},{"metadata":{"trusted":true},"cell_type":"code","source":"# compute power spectrum\naxis = 2 #we have placed the redshift-space distortions along the z-axis for the halos\nPk_h = PKL.Pk(delta_h, BoxSize, axis, MAS, threads, verbose)\n\n# Pk is a python class containing the 1D, 2D and 3D power spectra, that can be retrieved as\n\n# 1D P(k)\nk1D_h      = Pk_h.k1D\nPk1D_h     = Pk_h.Pk1D\nNmodes1D_h = Pk_h.Nmodes1D\n\n# 2D P(k)\nkpar_h     = Pk_h.kpar\nkper_h     = Pk_h.kper\nPk2D_h     = Pk_h.Pk2D\nNmodes2D_h = Pk_h.Nmodes2D\n\n# 3D P(k)\nk_h       = Pk_h.k3D\nPk0_h     = Pk_h.Pk[:,0] #monopole\nPk2_h     = Pk_h.Pk[:,1] #quadrupole\nPk4_h     = Pk_h.Pk[:,2] #hexadecapole\nPkphase_h = Pk_h.Pkphase #power spectrum of the phases\nNmodes_h  = Pk_h.Nmodes3D","execution_count":12,"outputs":[{"output_type":"stream","text":"\nComputing power spectrum of the field...\nTime to complete loop = 7.12\nTime taken = 12.30 seconds\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"Lets compare ths power spectrum with the one from matter:"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.xscale('log')\nplt.yscale('log')\nplt.xlabel(r'$k~[h{\\rm Mpc}^{-1}]$')\nplt.ylabel(r'$P(k)~[(h^{-1}{\\rm Mpc})^3]$')\nplt.plot(k, Pk0)\nplt.plot(k_h, Pk0_h)\nplt.plot([1e-1,3],[BoxSize**3/pos_h.shape[0], BoxSize**3/pos_h.shape[0]])\nplt.legend(['Matter', 'Halos', 'Expected shot-noise'])\nplt.show()","execution_count":13,"outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 640x480 with 1 Axes>","image/png":"iVBORw0KGgoAAAANSUhEUgAAAkIAAAG9CAYAAAD5ixlRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABon0lEQVR4nO3dd3hUZd7G8e+k94QQSIGE3ktCC00kaJSiCKjo2igq6ivYWAu4irpr2bViia5lARuCroiuHVGkiIBA6FVCQktCgPQ+M+8fB0ZiQkifSeb+XNe5kjlz5pzfxMHcec5TTFar1YqIiIiIE3KxdwEiIiIi9qIgJCIiIk5LQUhEREScloKQiIiIOC0FIREREXFaCkIiIiLitBSERERExGkpCImIiIjTcrN3AY7MYrFw9OhR/P39MZlM9i5HREREqsBqtZKTk0NERAQuLpW3+SgIVeLo0aNERkbauwwRERGpgUOHDtG6detKj1EQqoS/vz9g/CADAgLsXI2IiIhURXZ2NpGRkbbf45VREKrEmdthAQEBCkIiIiKNTFW6taiztIiIiDgtBSERERFxWgpCIiIi4rQUhERERMRpKQiJiIiI01IQEhEREafV5INQUlISI0aMoHv37vTq1Yu8vDx7lyQiIiIOosnPIzRlyhSefPJJhg0bxsmTJ/H09LR3SSIiIuIgmnQQ2rFjB+7u7gwbNgyA4OBgO1ckIiIijsShb42tXLmSsWPHEhERgclkYunSpeWOSUhIoG3btnh5eTFw4EDWr19ve27fvn34+fkxduxY+vbty9NPP92A1YuIiIijc+gglJeXR3R0NAkJCRU+v3jxYmbOnMljjz3Gpk2biI6OZuTIkaSnpwNQWlrKqlWreP3111m7di3Lli1j2bJlDfkWRERExIE5dBAaPXo0Tz75JBMmTKjw+RdffJFp06YxdepUunfvzr///W98fHyYN28eAK1ataJ///5ERkbi6enJmDFjSExMPOf1ioqKyM7OLrOJiIhI0+XQQagyxcXFbNy4kfj4eNs+FxcX4uPjWbt2LQADBgwgPT2dU6dOYbFYWLlyJd26dTvnOZ955hkCAwNtW2RkZL2/DxERkSattAh2fg47vwCr1d7VlNNog1BGRgZms5nQ0NAy+0NDQ0lNTQXAzc2Np59+mgsvvJDevXvTqVMnLr/88nOec/bs2WRlZdm2Q4cO1et7EBERafIKs+DjSfDxTVCF1eAbWpMeNQbG7bXRo0dX6VhPT08NrxcREalLDtgKdLZG2yIUEhKCq6sraWlpZfanpaURFhZWq3MnJCTQvXt3BgwYUKvziIiIyJkg5HitQdCIg5CHhwf9+vVj+fLltn0Wi4Xly5czePDgWp17+vTp7Ny5kw0bNtS2TBEREed2pkXIAW+LgYPfGsvNzWX//v22x0lJSSQmJhIcHExUVBQzZ85k8uTJ9O/fn9jYWObOnUteXh5Tp061Y9UiIiLyB8duEXLoIPTbb78xYsQI2+OZM2cCMHnyZBYsWMC1117L8ePHmTNnDqmpqcTExPDtt9+W60AtIiIiduLgLUImq9XBezHZUXZ2NoGBgWRlZREQEGDvckRERBqfrMPwUg9wcYc5GQ1yyer8/m60fYTqkzpLi4iI1BEHbxFSEKqAOkuLiIjUFcfuI6QgJCIiIvVHLUIiIiLivNQi1Oioj5CIiEgdU4tQ46E+QiIiInXEqhYhERERcVrqIyQiIiLOSi1CIiIi4vTUItR4qLO0iIhIHVGLUOOjztIiIiJ15UwfIftWcS4KQiIiIlJ/1CIkIiIizkujxkRERMRZqUVIREREnJdahERERMRZqUWo8dHweRERkbqiFqFGR8PnRURE6ohahERERMR5qUVIREREnJVahERERMR5qUVIREREnJVahERERMR5qUWo0dHweRERkTqiFqHGR8PnRURE6opahERERMRZnWkQUouQiIiIOJ8zLUL2reJcFIRERESk/qiPkIiIiDgv9RESERERZ6UWIREREXFeahESERERZ6UWIREREXFeahFqdDSztIiISB1Ri1Djo5mlRURE6opahERERMRZqUVIREREnJdahERERMRZFeUYX0uL7FvHOSgIiYiISP35eJLx9VSSfes4BwUhERERqT+WUntXUCkFIREREXFaCkIiIiLitBSERERExGkpCImIiIjTUhASERERp6UgJCIiIk5LQUhEREScloJQBep99fmCU7Drf3D4N8g6DOaS+rmOiIiIo2je0d4VVMhktdpWQ5M/yc7OJjAwkKysLAICAuruxMm/wPzRZ+0wgW8I+Ief3sIq/uobAi6udVeHiIhIfbJa4Ykg4/sbPoVO8Q1y2er8/nZrkIrkT0zQqj/kpEJuqjHrZt5xY0vdWsnLXMEv9NxB6cxXn2CHXdxORESciLn4j+8j6+kuSy0pCNlDm8EwbbnxvcUC+Scg55gRjM71NS8drGbIOWpslXH1AL+w08GoktDkFajAJCIi9ac474/v3X3tV0clFITszcUF/FoYW3jvcx9nPt1qdM7AdPr7/AwjgWelGFtl3LzP07oUBh5+4OYJ7t7g6mnUKyIiUhXFucZXNy9wdczI4ZhVSXmubhAQbmyVKS2G3LTKW5dyjkFhJpQWGKsBV2dFYFcPI0C5eRofbHevP763baeDU7n9Z44/85qzj/GEZm0gqI1aqUREmorEhcbX0kL71lEJBaGmxs0DgiKNrTIlBWVbkir6mptmNGtazX+8zlxsbEX1VH9Aa2g7FNpeAG2GQnB7BSMRkcZqxTP2ruC8FISclbs3BLcztvMxlxpp3rYVGUGqtOj047O+LznrmOrsLymAk79D9mHYutjYwLhN1+Z0MGp7gTH8UsFIRMTx5ab/8f2lT9qvjvNQEJLzc3UDVz/w9Kvf6xTnwaH1kLwGDq6BI78ZrVPb/2tsYIyaazPUaDVqcwG06KJgJCLiiM4eBT3kLvvVcR4KQuI4PHyhwwhjA6OV6PAGIxQdXG18n5sGO5YYG4BPyB+hqO1QaNFNHbpFRBzBgRXG114T7VrG+SgIieNy94Z2FxobGLfRjmw83WK0Cg5tMEbJ7fzc2AC8g6HNEOM2WuRAaNnNOI+IiDSs7Z8ZX8Oj7VvHeSgISePh7nW6I/VQGP6gMULu6CYjFB1cA4fWQcFJ2P2lsQGYXIwO1y27QcsextfQHsY+zdItIlJ/sg8bXx38j1EtsVGJeltiQ+qHuQSOJhrBKHkNHNlkBKOKuHlBSGcjFJ0JSaHdjc7Z5+tzZLUa0w/kpp/e0ow5nnLTjD5M/W82pgMQEXFWeRnwXAfj+5m7ICCiQS+vJTbEObm6G1O4Rw6AYTONwJKbDuk7jS3t9Nfju6Ek3+jI9+clTbyCoGV3IxQ1awuFWUbAsYWedGOW77Onjf+zxIUwcQE071CPb1ZExIH9/C/ja2ivBg9B1aUWoUqoRaiJslgg8+DpYLQL0ncY35/YX3bOpPPxCjRagHxbgl9L8GkO2z81WqE8/OCyFyD6L/X2NkREHNbc3pCZDNctgi6jz398HVOLkEhlXE73GwpuD90u/2N/aRFk7DXCUdoOyDpkdL72Ox10zg49vi2MPkt/NmwmLLnNuD332e3GqIkxz9f/1AMiIo7CajUm5QVo0dW+tVSBgpDIGW6eENbL2GoqIAImfQ6rXjBmVN3ykTE30tXzICKmzkoVEXFYB1aAuchYRsnBb4sBaMIVkbrm4mqMapvyFQS0MmbM/s8l8Osbxl9KIiJN2f4fjK+dL20UA0eafBBq27YtvXv3JiYmhhEjRti7HHEmbYbAHauh6+VG5+pvZ8FH10HeCXtXJiJSf3Z/ZXztPs6+dVRRkw9CAL/88guJiYn89NNP9i5FnI1PMFz7gdFPyNUT9n4D/x5qzJQtItLUHNkIp5LA5Aod4+1dTZU4RRASsSuTCWKnwbTl0LyTsX7au2Phx6egON/e1YmI1J2fnja+th1qjKxtBBw6CK1cuZKxY8cSERGByWRi6dKl5Y5JSEigbdu2eHl5MXDgQNavX1/meZPJxPDhwxkwYAAffvhhA1UuUoGwXnD7z9DnRrBaYOWz8FIP+PFJyEmzd3UiIrVjtRojbgE6j7JvLdXg0EEoLy+P6OhoEhISKnx+8eLFzJw5k8cee4xNmzYRHR3NyJEjSU9Ptx2zevVqNm7cyBdffMHTTz/N1q1bKzyXSIPw8IVxCXDVfyCojTHn0MrnYG5PWDr9j/+JiIg0NqnbjBZvdx/of4u9q6myRjOhoslk4rPPPmP8+PG2fQMHDmTAgAG89tprAFgsFiIjI7nrrruYNWtWuXM88MAD9OjRgylTplR4jaKiIoqKimyPs7OziYyM1ISKUj8sZqNT4drXjHXSzmg/AgbPgI4Xn3+5j+J8Y9Iy3xbgG1K/9YqIVGbZY7BmLnQeDdcvsmspTjGhYnFxMRs3bmT27Nm2fS4uLsTHx7N27VrAaFGyWCz4+/uTm5vLjz/+yDXXXHPOcz7zzDM88cQT9V67CGAMs+9+hbEd2gC/JsDOz+HAT8bWoisMnm40MWcdgpNJxnbq9NeTByD39KRlJlcjOPW+Frpe5vCLHIpIE2OxQOLp7idnT1TbCDTaIJSRkYHZbCY0NLTM/tDQUHbv3g1AWloaEyZMAMBsNjNt2jQGDBhwznPOnj2bmTNn2h6faRESqXeRAyByAZxKhnVvwqb3jDXRvrjr/K/18IfiHNj3vbF5BhjhKvo6iBpizKQtIlKfklYYi0+7+0Kvczc4OKJGG4Sqon379mzZsqXKx3t6euLp6fiTP0kT1qwNjHoa4h4ywtCv/4bsw+AXBsHtjGVBmrU7/X0743ufYMjYD1sXwdbFkJkCmz8wtsAo6H0N9JsMQVH2fnci0lStecX42ucGcPOwby3V1GiDUEhICK6urqSllR1tk5aWRlhYmJ2qEqkjXoEw5C6jr1BpUcXrmp0tpCNc9AjEPQwpa41QtGMpZKXAqudhzcvQbwpceD/469+HiNShwuw/+jn2nWTfWmqg0baZe3h40K9fP5YvX27bZ7FYWL58OYMHD67VuRMSEujevXult9FEGoTJdP4QdDYXF2P+jitehfv3GmuctR0GlhLY8Da8HA3fP6LZrUWk7qx5GUryjXnSWvawdzXV5tCjxnJzc9m/fz8Affr04cUXX2TEiBEEBwcTFRXF4sWLmTx5Mm+++SaxsbHMnTuXjz/+mN27d5frO1QT1el1LuLQklbC8n/A4dPzbHn4waA7jc7Y3kF2LU1EGrHsY/Byb2MZoWveN/onOoDq/P526CC0YsWKCtcHmzx5MgsWLADgtdde47nnniM1NZWYmBheeeUVBg4cWCfXVxCSJsVqhX3L4Md/QOrp+bS8gqDLGAjrCaGnN9/mdi1TRBqRVS/A8r9Dq/5w6w/nn/KjgTSZIGRvCkLSJFmtsOsLYyr847vLP+8fbgSiVv2MWbCDNHJSRCpQmAVze0NhJlz2Agy41d4V2SgI1VJCQgIJCQmYzWb27t2rICRNk8UM+5fD0U2Qth1StxtzFJ3N5GrMCTLw/yBqkMP8tSciDuCr+42+h/7hcM8WcHOcUdcKQnVELULidIpyIG0npG0zJndMWvnHc+ExMOj/oMcEh/ofnojYQfJamH96PbGr/gO9rrZvPX+iIFRHFITE6aXtgHX/hq0fQ2mhsc8v1BjW3/9m8PSzb30i0vDMpfCP030J/cPh3u3g6liz8VTn93ejHT4vIg0gtIcxFP++nXDRo8b/9HLTYNmjxkiRVS8Yc4iIiPPYetY6YuMSHC4EVZeCUAU0j5DIn/g2NyZjvHcbjHvdmOE6/4QxWmRuL1jxLyjItHeVIlLfUrfD59ON73tdY6xx2Mjp1lgldGtM5BzMpbD9U2PW6oy9xj4Pf+g/1ZifKCDcvvWJSN3LSYX3xsPxXUbr8N2bHXaBZ/URqiMKQiLnYTHDzqWw8nlI32nsc3GH6GthyD3QorNdyxOROmAxw+b3jZbfnKOACab9CK362ruyc1IQqiMKQiJVZLXCvu9h9VxI+eX0ThN0GwsXz4GQTvasTkRqylwCi2+Evd8aj/3CYNJSaNnNrmWdj4JQHVEQEqmBQ+uNtYd2f2k8NrkaCzHGzdKCryKNibkEFt8Ee78xHg/8PxgyAwJb27euKlAQqiMKQiK1kL7b6Ey95yvjsbuPsbbZkLvBS/+eRBzempdh2RzABNctgi6j7F1RlWn4fC1p1JhIHWjZFa5bCFO/hdYDjNWpVz4HL/WE7/4GmSn2rlBEzuXQemOhZoBhf21UIai61CJUCbUIidQRqxV2/c9Y8PXMKDOTq9GHKHaa0dTu5g3uXkbLkau7fesVcWaF2fDvCyAzGbqPg4nvNrrldXRrrI4oCInUMYsF9i+DtQmQ9PO5jwuPgd7XQM+r1K9IpCFZzPCfS+HIb+Dbwhgi7+lv76qqTUGojigIidSj1O3w6xuw7zsozjNunf2ZyQXaXQhth0FQFARGQlAkBLRqdH+hijQKq16E5U8YLbM3fwvh0fauqEYUhOqIgpBIA7JaobQICk7C7q9g2ydwaF3Fx7bsAcNmQvfxjX56fxGHsX85fDgRrGYY/SwMvN3eFdWYglAdURASsbNTB2HHUji+B7IOnd4Og6XUeL5ZOxhyF7QZYiz74eZpz2pFGq/U7fBOPJQWQK+JcOXbjbrVVUGolhISEkhISMBsNrN3714FIRFHUpAJG96Gta8brUdnmFygWVujc+dFj4KLq70qFGlcTiXDvJGQcwwiB8JNn4GHr72rqhUFoTqiFiERB1acBxvfhe3/hYx9UJT9x3Pdxxt/0bp52K08kUYhJw0WjIET+yEwCm5bYSyy3MgpCNURBSGRRsJqhdw02P8D/O9esJRAx0vgmvfAw8fe1Yk4prSd8OaFxr+XwEiY+o0xGKEJqM7vb/UyFJHGz2Qyhtn3udFYFXvRDcYw/ffGGf0dmrc3+hO5uhuhyeRijDxz0Zyy4qQOroaPrjdCkFcg3LikyYSg6lIQEpGmpePFRh+HhdfA4fXGVpFm7aD/VIi5AXxDGrZGEXuxWGDze/D1A2AuNvoEXT2vUawfVl90a6wSujUm0ogd3wub34cTv8PJ340lPawWwGT8FXxm5JmLO3g3M1qJXD0gcgB0HgUd48En2K5vQaROnToIn8+Ag6uMx10vh6veAXdvu5ZVH9RHqI4oCIk0UcV5sP1T+G0eHN187uNadIXIWGgda/zl3LyjbqdJ41NaBBvegR+fgpI8I/APuRtGPNxkR1cqCNURBSERJ3AyyQhGWI2h+Qd+gr3fQdr28sd6Bf0RjNrHQev+jXquFWniSoth1xfw45NwKsnY12YojHvNmHerCVMQqiXNIyQi5GUYK3AfXm98PbIRSgvLHtOyO8Rcb3Q2NRcbv1yihhiLx4rYS04qbFxgtHjmphn7/MJgxGzoM8kpWjUVhOqIWoRExKa0GNK2waENkLLWaDUqLSh/nLuP0WF7yN1G65FIQyjMhj3fGLd8f1/+Rx84v1DofzMMnt4oF0+tKQWhOqIgJCLnVJAJWz82huljMvpaHNkEual/HNNuOIx6BkJ7GI+tVijMMmbEDmilJUGkdorzjQEBe76B5DVGq+QZkYMgdhp0u8IpJxZVEKojCkIiUi1WK6Rug/VvwpZFxl/lJlfoexNkHzN+WRXnGse6uBm31gbPgN7XqK+RVE9mCiy63vi8neEXBp0vhQG3NtpV4+tKvQehL774otpFXXLJJXh7N64hegpCIlJjmSnw3d+Mzqp/5uoJ5qI/HncfbwzZ9wk2bqd5N2uwMqUR2v4pfHX/6bX2TBB7mxG2w3rZuzKHUe9ByKWaHa1MJhP79u2jffvG1UtdQUhEam3PN7BjKbToYsxN1LyjMW9L1mFIXAg//wus5j+ON7kaYahjvDFk36e5sZislgqRvBPw/d9gy0fG4/BouPZDp50RujINEoRSU1Np2bJllY739/dny5YtCkIiIn92eCOs+zfkZxjhKGNv+WNc3KBFN+NWm9UCYT2hzRBjQjz/sIavWRqWxQJbF8H3j0D+CWytQJc80SQnQ6wL9b7W2OTJk6t1m+vGG29UkBARqUjrftD67T8eZ6bAvmXGdmKfMYy/MNMYsXZGxp4/bo9EDjSWRzixHzpdAhfMVOtRU3J4I3x5zx99gVr2gCteMeawkjqhztKVUIuQiNid1Wq0FKVtBzcvo1XoaCLs/RaO/Fb+eA9/8A81Jn8MioT2I4xRa4fWGaPaWnaDLmMgtHtDvxOpKosZ9i+HXxPgwApjn2cgDLvP6Fzv6m7X8hqDer01durUKaxWK8HBwRw/fpxVq1bRpUsXevToUauiHZGCkIg4tKzDsPtrY0i+b3NY+QJkH67aa4OijL5HFz1qLLlQUmB00nZxNfYX5RjBS61LDcNcaowq3Put0acs5+gfz/W+FkY+rcWBq6HegtA777zD008/DcADDzzAhx9+SHR0NCtXruSee+7h1ltvrV3lDqK+Z5YuLDGTcjKfqGAfvNyb5jovImIHpcXGLbLCTMg/adxOObja2OcfZoxMO7bFmPvozIR7FTIBp381hPWGfpONodnezYxfxs07GcP9NeS/dvJOwLHNsOm90xN0njVzuVcg9LnJ6AvUrI39amyk6i0I9e7dm3Xr1lFQUEBUVBRJSUm0aNGCrKwshg8fTmJiYm1rdyj11SK0OeUUE17/BZMJIgK9adPch7YhvrRr7kub5j60C/ElUiFJROpL3gmj/9H6t4y+Ru6+RsApyjbCVEle5a/3DDBakMJ6QWArCIwyOm9HDoSsQ8btPP8w8Gup2zhnHN8LlhLI2Gcs13JwVfkFf919oecE6DQSOl2qpVpqod46S7u5ueHt7Y23tzcdO3akRYsWAAQGBmLSXwZVlplfgr+nGzlFpRzJLOBIZgG//H6izDFnQlLbEB/aNvc1thBf2oX4EBnsg6ebQpKI1JBvc2OLHAjDZxmdrc/cArNYjNmxvYONYLRxgdG/qCDTaGnKPmbsBzi6ydjA6M/yZyYXiOhrLO1QcMoY7ebbEjz9jBXRvYKMfbnpxkzb3s2M6QL8WtT/z6CuWSxQnAMnfjfCzol9Rotc2g5j1N+ZRU//zCvQGP3X9XLoPLLJrgbvyKrVIjRgwABWrVqFl5cXWVlZBAYGApCbm8uwYcPYvHnzec7QuNRnHyGr1crJvGIOnsgjKSOf5BN5JGXkcfBEHgcz8sktOnez9ZmQ1C7kjxYkIygpJIljsVqtFFS0Hpc0XpZSSN9tLA+SstZY1uHEfkj51Rj67+FndOzNSzvP7bdzMLkanbwLs4xFbE0m8AwyAoKnvzHnkn+E8dUn2JhO4NRB8A011n4L6WRc3yvAWIKiKNsIXCGdjL5QJhfw8DbqzDxkDD/PPmoc5+JmnNdigdxjUFxgvM+Ck1CQBViNMGgFUrdC/inwbwm5xyEvvWrvN7QntOoLETHQ9kIIiKj+z6gJ8nbzrtMGlXq7NXbmhH8uNj09neTkZAYMGFCzih2UvTpLW61WTuQVczDDCEfJJ/JJOpHHwQxjyys2n/O1LiaICPK2BaOzW5Oign3wcGv6qw6L48gvyWfgwoH2LkNEHNy669fh4153HfPr7dbYmRagP2vZsmWVJ1eU8zOZTIT4eRLi50n/tsFlnrNarWTkFp9uOfqjBenM47xiM4dPFXD4VAGr95c975mQdKYFydaaFOJLZDOFJBERcT51Mo9QSUkJqamp5Ofn06JFC4KDg8//okagsQ2ft1qtHM8tMlqQTrce2b4/kUf+eVqSWjXzLtOC1OZ0h23r6dEjZz4pZ39gznx8bPtsx/zpNX96bbnXAW2a+9Al1F/9zZoQ3RoTuyktNNZ0O/P/k5J8ONPikJtudPYOjARLsTFNgNhVo7k1dracnBw++OADFi1axPr16ykuLsZqtWIymWjdujWXXnopt912W6O+XdbYglBlzoSkgxn5f7QkndU/qbKQ1JDahfgyumcYY3qF0yOi/G1YERGR86n3IPTiiy/y1FNP0aFDB8aOHUtsbCwRERF4e3tz8uRJtm/fzqpVq1i6dCkDBw7k1VdfpVOnTjV+Q/bSlIJQZaxWK8dzisr1Rzp0Kp9Ss/HxOBNIzsSSM/nk7JxiOv2s7TnKHlzutbbHJswWKzuPZVNcarGdLyrYh9G9whjTM5zerTUyUUREqqbeg9B1113HI488ct7ZpIuKipg/fz4eHh7cfPPN1b2M3TlLEHIUuUWl/Lg7nW+2HeOnPekUlvwRiloFeTOmVxije4XTJzJIoUhERM6pQW6NOQMFIfvJLy7lp93H+Xr7MX7clU5ByR+37iICvRjVM5wxvcLoG9UMFxeFIhER+YOCUB1REHIMBcVmft6bztfbUlm+K63M9AGhAZ6M7hnO6J5h9G8bjKtCkYiI02vQIPTMM88QGhpa7tbXvHnzOH78OA899FBtTm9XCkKOp7DEzMq9x/lmeyo/7Ewj56yJJ1v4ezKqRxije4UR2zYYN1dNByAi4owaNAi1bduWhQsXMmTIkDL7161bx1/+8heSks4xrbgDq+9FV6VuFJWaWb0vg6+3pbJsZyrZhX+Eoua+HlzUtSXRkUH0bh1IlzB/zbgtIuIkGjQIeXl5sWvXLtq1a1dm/4EDB+jevTuFhYXneKXjU4tQ41FcamHN7xl8s+0Y3+9MIzO/pMzz7q4muoYF0Kt1IL1bBdKrdSCdQ/1xV6uRiEiTU28zS1ckMjKSNWvWlAtCa9asISJCa6hIw/Bwc2FEl5aM6NKSp8wW1v5+gnVJJ9h6OIttR7LIzC9h2xHj+4VnvaZ7eAC9WwfSq1UgvVsH0aGFb5VuqVmtVrILSjmeW8SJ3CJO5BWTkVtES39PLu0epg7cIiKNRK2D0LRp07j33nspKSnhoosuAmD58uU8+OCD/PWvf611gSLV5e7qwoWdW3BhZ2MFa6vVyuFTBWw9nMXWI5lsOx2OcgpLSTyUSeKhTNtrvd1d6RFhtBx1bOlHbmGpLeRk5BYboSe3mBN5RZSYK25MvbhrS56fGE0zX4+GeLsiIlILtb41ZrVamTVrFq+88optdmlvb28eeugh5syZU1d12oVujTVdFouV5JP5bD1sBKOtR7LYcSSr0gVtK+Lv5XZ6XTgPAr09WLnvOMWlFiICvXj1+r70a9Osnt6BiIici12Gz+fm5rJz5068vb3p3Lkznp6edXFau1IQci5mi5WkjFyj5ehwFskn8gjy8aC5rwfNT4edM4vhNvfzINjXAy/3sh2wdxzNYvqHmzh4Ih83FxMPjurCrRe0160yEZEG1OBB6D//+Q8vvfQS+/btA6BTp07ce++93HrrrbU9tV0pCElN5BSWMHvJNr7cegzQrTIRkYZWnd/ftR4yM2fOHO655x7Gjh3LJ598wieffMLYsWO57777Gv2tMZGa8Pdy59Xr+vDUhJ54uLmwfHc6l72yio3Jp+xdmoiI/EmtW4RatGjBK6+8wnXXXVdm/0cffcRdd91FRkZGrQq0J7UISW3tOJrFjIWbScrIw83FxAMjuzBtmG6ViYjUpwZtESopKaF///7l9vfr14/S0tIKXiHiPHpEBPLFjKGMjY6g1GLlmW92c+t7v3Eqr9jepYmICHUQhG666SbeeOONcvvfeustbrjhhtqeXqTR8/dy55W/xPD0hF54uLnw4+50xryyio3JJ+1dmoiI06v1rbG77rqL9957j8jISAYNGgQYy2ukpKQwadIk3N3dbce++OKLtau2genWmNS1nUezmb5wE0kZebi6mLj1gnbcE98JH49aT+klIiKnNeiosREjRlTpOJPJxI8//libSzU4BSGpD7lFpfzts218nngUgFZB3vx9XA8u7hZq58pERJoGu8wj1BQpCEl9Wr4rjTmf7+BIZgEAo3qE8dgV3QkP9LZzZSIijVuDdpYWkZq5uFsoy2ZeyO0XtsfVxcS3O1KJf+Fn5q1OotRssXd5IiJOocYtQjfffHOVjps3b15NTu8Q1CIkDWXXsWz+9tk2NqVkAtCzVQBPje9FdGSQXesSEWmMGuTWmIuLC23atKFPnz5UdorPPvusJqd3CApC0pAsFiuLNhzin9/sIruwFJMJJg1qw32XdMbN1YWiEjOFpRYKS8ynNwtFJWaKSi10aOFHVHMfe78FERGH0CBBaPr06Xz00Ue0adOGqVOncuONNxIcHFyjgh2VgpDYw/GcIp7+ehefbT5S5deYTHBJt1BuH96efm2a1r9DEZHqarDO0kVFRSxZsoR58+bxyy+/cNlll3HLLbdw6aWXYjI5zsy5+fn5dOvWjYkTJ/L8889X+XUKQmJPa/Zn8OjS7RzIyLPt83B1wdPdBS93V7zcXfByc8XFZGJPWo7tmP5tmnHbhe2J7xaqGaxFxCnZZdRYcnIyCxYs4L333qO0tJQdO3bg5+dXF6eutb/97W/s37+fyMhIBSFpVKxWK1kFJXi6ueLp5nLOYLMvLYe3Vh5gaeIRSszGP+n2LXyZNqw9E/q0wsvdtSHLFhGxK7uMGnNxccFkMmG1WjGbzXV12lrbt28fu3fvZvTo0fYuRaTaTCYTQT4eeHu4Vtq60ynUn+cmRrP6oYu4Y3gH/L3cOHA8j9lLtnHBv37ig1+TK+3LJyLirGoVhIqKivjoo4+45JJL6Ny5M9u2beO1114jJSWlTlqDVq5cydixY4mIiMBkMrF06dJyxyQkJNC2bVu8vLwYOHAg69evL/P8/fffzzPPPFPrWkQag9AAL2aN7sra2RfzyGXdiAj0IiO3iEeWbufmBRs4nlNk7xJFRBxKjYPQnXfeSXh4OP/85z+5/PLLOXToEJ988gljxozBxaVuGpry8vKIjo4mISGhwucXL17MzJkzeeyxx9i0aRPR0dGMHDmS9PR0AD7//HM6d+5M586d66QekcbCz9ONW4e15+cHR/Do5d3xcHPhpz3HGTV3JT/sTLN3eSIiDqNWw+ejoqLo06dPpR2jlyxZUuPizmYymfjss88YP368bd/AgQMZMGAAr732GgAWi4XIyEjuuusuZs2axezZs/nggw9wdXUlNzeXkpIS/vrXvzJnzpwKr1FUVERR0R9/MWdnZxMZGak+QtLo7UnN4Z5Fm9mdanSqvn5gFI9c1k1rnIlIk9QgfYQmTZrEiBEjCAoKIjAw8JxbfSkuLmbjxo3Ex8fb9rm4uBAfH8/atWsBeOaZZzh06BAHDx7k+eefZ9q0aecMQWeOP7v2yMjIeqtfpCF1CfPn8xlDmTasHQAL16Vw+Sur2Xo4076FiYjYWY3/HFywYEEdllF9GRkZmM1mQkPLLlQZGhrK7t27a3TO2bNnM3PmTNvjMy1CIk2Bp5srf7usO3FdWvLXj7dwICOPK1//hcEdmhMR6E14kJfta3igN62beWu0mYg0eU7TLj5lypTzHuPp6Ymnp2f9FyNiR0M7hvDtvcP429LtfLX1GKv2ZVR4nJ+nG09c0YOr+rVu4ApFRBpOnQWhdevWMXDgwLo63XmFhITg6upKWlrZjp9paWmEhYU1WB0ijVGQjwevXdeHWy9ox/70XI5lFXIsq4CjmX98zS0q5a+fbGHDwZM8fkUPtQ6JSJNUZ0Fo4sSJpKSk1NXpzsvDw4N+/fqxfPlyWwdqi8XC8uXLmTFjRq3OnZCQQEJCgkPNhyRS10wmE32imtEnqlm558wWK6/+uI+Xl+9j0YZDbDmcxes39KVdiK8dKhURqT/VGjV2zTXXVLjfarXyzTffkJubW2eFAeTm5rJ//34A+vTpw4svvsiIESMIDg4mKiqKxYsXM3nyZN58801iY2OZO3cuH3/8Mbt37y7Xd6gmNLO0OLvV+zK4Z9FmTuQV4+fpxnNX92Z0r3B7lyUiUql6W2IjODiY999/v9xkiVarlWuvvbbcbaraWrFiBSNGjCi3f/LkybbO2q+99hrPPfccqampxMTE8Morr9TZLToFIRFIzSrkro82seHgKQCmDm3Lw2O64e5aZxPTi4jUqXoLQldeeSX33nsvF154YbnnLrnkEpYtW1b9ah2YgpCIodRs4fnv9/Lvn38HYEiH5rxxQz8CfdztXJmISHl2WXS1KTm7j9DevXsVhEROW7YzjXsXbSav2EyHFr7MnxJLVHMfe5clIlJGgwWh1NTUJj1CSy1CIuXtPJrNLe9u4FhWIcG+Hrx1Uz/6tw22d1kiIjYNtvr8pZdeWpuXi0gj1D0igM+nD6VXq0BO5hVz/dvr+DzxiL3LEhGpkVoFId1VE3FOLQO8WHz7IC7tHkqx2cI9ixJ54fs9lJot9i5NRKRaahWEKltsVUSaNh8PN/59Yz9uu7A9AK/+uJ9r3/qVlBP5dq5MRKTqNP61AgkJCXTv3p0BAwbYuxQRh+biYuLhMd14+S8x+Hu6sTH5FGNeWcV/Nx5Wi7GINAq16izdu3dvtm7dWpf1OBR1lhapukMn8/nrx1tYf/AkAJf1CuepCT0J8vGwc2Ui4mwarLO0q6vWHhIRQ2SwDx/dNogHRnbBzcXEV9uOcdkrq9mTmmPv0kREzqlWQWjz5s11VYeINAGuLiamj+jIkjuH0La5D0cyC7j6jV9Ys7/iFe5FROxNfYREpM71bh3E0ulDiW0XTE5RKZPnree/Gw/buywRkXIUhCqgztIitRfk48H7t8RyRXQEpRYr93+yhbk/7FUnahFxKLVeYmPDhg3MmjWL48eP07FjR2JiYmxbVFRUXdVpF+osLVJ7FouV57/fw+srjHXKLu8dzt/H9STYV52oRaR+NOhaY127diUqKoorrriCpKQkEhMTSUxM5NSpUzRr1owTJ07U5vR2pSAkUncWrkvh0c+3Y7ZYCfb14LGx3bkiOkLzkYlInavO72+32l7s0KFDfPXVV3To0KHM/uTkZBITE2t7ehFpIq4fGEX3iABmfbqV3ak53LMokc82H+HJ8T1p3UwLt4qIfdS6RejSSy/lb3/7G8OHD6+rmhyGWoRE6l5xqYW3Vv7OK8v3U2y24OXuQnTrILpHBNAtPIDu4cZXVxe1FIlIzdT7rbErr7yS3r17Ex0djdVq5fXXX+eTTz6hWbNmNS7aESkIidSf/em5PLxkm20CxrO1CvJmypC2XBsbSYCXux2qE5HGrN6D0AMPPEBiYiJbtmwhI8OYH6R58+aMGzeOQYMG0adPH3r16oWHR+PsDJmQkEBCQgJms5m9e/cqCInUE6vVyu7UHHYczWbXsWx2Hs1m+5EscopKAfD1cGVi/0imDm1Lm+a+dq5WRBqLBu0sfeTIEVsH6TPbgQMHcHNzo0uXLo16CQ61CIk0vMISM0s3H2HemiT2puUCYDLBJd1CueWCdsS2C1YHaxGpVIMGoYrk5ubaWoymT59e16dvMApCIvZjtVpZtS+D/6xO4ue9x237e7YK4P5LuxDXpaUdqxMRR1bvQSglJaVacwQdOXKEVq1aVfcydqcgJOIY9qfnMG/NQZZsOkxhiQWA/4vrwF8v6Yybq+aFFZGy6n3R1QEDBnD77bezYcOGcx6TlZXF22+/Tc+ePfn0009rchkREQA6tvTn6Qm9WDvrYiYNbgPAGyt+5/q315GWXWjn6kSkMatRi9CJEyd46qmnmDdvHl5eXvTr14+IiAi8vLw4deoUO3fuZMeOHfTt25dHH32UMWPG1Eft9U4tQiKO6autx3jo063kFpXS3NeDW4e156p+rWjp72Xv0kTEATRYH6GCggK++uorVq9eTXJyMgUFBYSEhNCnTx9GjhxJz549a3pqh6AgJOK4DhzP5c4PN7E7NQcwVr4f0aUlt13Ynth2wXauTkTsye6dpZsKBSERx1ZYYubzxCMs3nCITSmZgBGIXro2hiuiI+xbnIjYTb33EWrqtPq8SOPg5e7KtQOiWHLnUJbddyGX9Q7HbLFyz6LNfPzbIXuXJyKNQI1ahLZu3UrPnj1xcalajtqxYwddunTBza3WS5s1KLUIiTQuFouVRz/fzofrUgD4x7ge3DS4rX2LEpEGV+8tQn369KnWqvKDBw8mJSWlJpcSEakyFxcTT47vyc1D2wHw6Oc7+Mtba/lpdzrqBSAiFalRE43VauXRRx/Fx6dqK0YXFxfX5DIiItVmMpl49PJuBHi78dqP+/n1wEl+PXCSLqH+/O2yblzYuYW9SxQRB1KjW2NxcXHVnuJ+4cKFhIeHV/dSdqVbYyKN29HMAuavSWLhuhTyis0A3DgoiofHdMPHo3HdqheRqtOosTqiICTSNGQVlPDi93t4d20yAG2a+/DqdX3o3TrIvoWJSL3QqDERkbMEervzxLiefHjrQCICvUg+kc91b/3KugNV7+soIk2TgpCIOI2hHUP49r4LGdKhOXnFZibPX8/qfRn2LktE7EhBSEScSoCXO/OmDCCuSwsKSyzc/O4GXvtxHwcz8uxdmojYgfoIVUJ9hESarqJSM3d/tJnvdqTZ9nUPD2DK0LZc2aeVVrUXacTUWbqWEhISSEhIwGw2s3fvXgUhkSaq1Gzh002H+XLrMX75/QRmi/G/w3YhvtxzcSfGRkfg6lK9EbIiYn8NHoRKSkpITU0lPz+fFi1aEBzcNBY8VIuQiPM4lVfMx78d4s2VBziZZ8x91rGlH/fGd2JMz3BcFIhEGo0GCUI5OTl88MEHLFq0iPXr11NcXIzVasVkMtG6dWsuvfRSbrvttka9XpeCkIjzyS0q5d1fDvLWygNkFZQA0LNVAM9M6E2v1oF2rk5EqqLeg9CLL77IU089RYcOHRg7diyxsbFERETg7e3NyZMn2b59O6tWrWLp0qUMHDiQV199lU6dOtX4DdmLgpCI88ouLGHe6iT+syqJnKJSXExwywXtuO+SzpqMUcTB1XsQuu6663jkkUfo0aNHpccVFRUxf/58PDw8uPnmm6t7GbtTEBKR4zlF/P3Lnfxvy1EAWvp78sDILlzVt7Vul4k4qAbtI5Senk7Lli1rcwqHpSAkImf8uDuNx7/YScrJfAA6h/px7YAoxsdE0NzP087VicjZGjQIXXjhhfz000+4urqWe660tBQ3t8bbhKwgJCJnKyo1s2DNQV77cT85RaUAuLuamHlJF+4Y3r7aazCKSP1o0CU2goKCuPvuu8vtP3HiBPHx8bU9vYiIw/B0c+X24R1Y/dBF/GNcD3q3DqTEbOVf3+5m+sJN5J0ORyLSeNQ6CL333nssW7aMefPm2fbt2rWL2NhYfH19a3t6ERGHE+jjzk2D2/L59KE8NaEn7q4mvt6WyvDnfuLZb3dzJLPA3iWKSBXVyTxC27ZtIy4ujm+++YZTp05x7bXXcsstt/Dcc8/h4tJ4Z2fVrTERqYqNyae4+6PNtgDk5e7CGzf0Y0TXptl/UsTR1XsfoSuvvJKYmBjbFhUVxUcffcRdd91FYWEhr776KlOnTq3xG3AUCkIiUlUlZgs/7Ezj7VUH2JSSiZuLiZeujWFsdIS9SxNxOvXeR6hDhw6sWrWKW2+9lbZt29K8eXPefvttrFYr119/PX379qWkpKRGxYuINEburi6M7hXO4tsHc0V0BKUWK3cv2sw9izbz+/Fce5cnIudQ61tjR44cITExscx24MAB3Nzc6Nq1K1u2bKmrWhucWoREpCbMFiv/+HInC345CICLCaYMaccDI7vg7VF+hK2I1C27L7qam5tLYmIiW7ZsYfr06XV9+gajICQitbH9SBYvL9/Hsp3GCvftQny5tHsoEUHejO4ZRssALztXKNI01XsQSklJISoqqsrHHzlyhFatWlX3Mnaj1edFpC79tDudWUu2kpZdZNvn6+HKPfGdmDq0He6ujXdQiYgjqvcgFBoayvjx47n11lvPuahqVlYWH3/8MS+//DK33XZbhXMNOTq1CIlIXcnKL+GzzYc5dKqA9Ukn2XYkC4ARXVrw5k398XBTGBKpK/UehE6cOMFTTz3FvHnz8PLyol+/fkRERODl5cWpU6fYuXMnO3bsoG/fvjz66KOMGTOmxm/GnhSERKQ+WCxW/rvpMHM+305hiYXRPcN4ekIvmvl62Ls0kSahwfoIFRQU8PXXX7Nq1SqSk5MpKCggJCSEPn36MHLkSHr27FnTUzsEBSERqU8/7z3OtHd/o9hsASDEz4ObL2jHzUPb4eWuTtUiNdVgQSg5OZmtW7cSGhpKbGxsTU/jsBSERKS+/bAzjSe+3MGhk3/MRt0qyJtZo7tyee9wrV8mUgMNEoQ++ugjpkyZQklJCSaTiT59+vDNN9/QokWLGhXtiBSERKSh5BeX8u32VJ79dg+p2YUA9IgIYHxMK8b3aUULf61wL1JVDbLo6hNPPMH111/P7t27+f777wGYNWtWTU8nIuLUfDzcuLJva368fzj3xXfG292VHUezeerrXYyau5LEQ5n2LlGkSapxi5CHhwd79+6lbdu2AOzevZt+/fqRl5dXl/XZlVqERMRejucU8c32Y3zwazJ703Lxcnfhsl4RxLZrxtX9InF10S0zkXNpkBah0tJSfHx8bI+7du2KxWIhNTW1pqcUEZHTWvh7MmlwWz67cyjDO7egsMTCp5sO89Cn2/j7/3ZQD3PhijilWk1c8e677/LLL7+Qm2uso+Pm5kZ+fn6dFCYiIuDr6cZ/JvfnP5P7c/vw9phM8O7aZP7x5S6OZhac/wQiUqka3xobPnw4iYmJ5OTk4OLiQrt27Th48CAPPvgg8fHx9O/fH39//7qut0Hp1piIOJp3Vh3gya922R5HBnszpH0I98R3IiLI246ViTiOBl1rbN++fWzcuJFNmzbZtszMTFxcXOjUqRO7du06/0kclIKQiDii/205yge/JrMu6aRtn5e7Cw+P6cZNg9poyL04PbsvupqUlMRvv/3G5s2befrpp+v69A1GQUhEHFlmfjFbD2fx2o/7WX/QCEVX9mnFzRe0o0dEgAKROC27B6GmQkFIRBoDq9XKWysP8M9vd3Pm/+i9WgUy46KOXNItFBeNMBMnoyBURxSERKQxWfv7Cd5be5Cf9qRTWGIs29E1zJ8ZF3VkdM9wDbkXp6EgVEcUhESkMTqRW8R/Vifx3tpkcotKAejY0o9Jg9vg6ebCiC4taRngZecqReqPglAdURASkcYsK7+EeWuSmL8miezCUtv+lv6eLL59MO1CfO1YnUj9aZAJFRuDzMxM+vfvT0xMDD179uTtt9+2d0kiIg0m0Med+y7pzOpZF/HAyC4M7dic1s28Sc8pYuK/f+G/Gw9jsehvYXFuTbpFyGw2U1RUhI+PD3l5efTs2ZPffvuN5s2bV+n1ahESkabmeE4RN76zjj1pOQBERwbxxBU9iIkMsm9hInVILUKnubq62pYBKSoqwmq1alp6EXFqLfw9+eKuocwa3RVfD1e2HMpkwutruG9xIusOnKCwxGzvEqUJOnQy32E/Ww4dhFauXMnYsWOJiIjAZDKxdOnScsckJCTQtm1bvLy8GDhwIOvXry/zfGZmJtHR0bRu3ZoHHniAkJCQBqpeRMQxebq5csfwDvx0fxwT+rTCaoXPNh/h2rd+JfqJ73nx+z0UFDvmLy1pfLYdzmLYsz8x5uVV9i6lQg4dhPLy8oiOjiYhIaHC5xcvXszMmTN57LHH2LRpE9HR0YwcOZL09HTbMUFBQWzZsoWkpCQWLlxIWlpaQ5UvIuLQWgZ48dK1MXw+fShX9m1FiJ8HRaUWXvlxP0P+uZw5n28n8VCmvcuURu6LLUcAOJCRZ+dKKtZo+giZTCY+++wzxo8fb9s3cOBABgwYwGuvvQaAxWIhMjKSu+66i1mzZpU7x5133slFF13E1VdfXaVrqo+QiDgTq9XKN9tTeeabXRw6aSzo6mKCey7uzNCOzekb1UyTM0q1Pf7FDhb8chCAg/+8rEGu6RR9hIqLi9m4cSPx8fG2fS4uLsTHx7N27VoA0tLSyMkxOgRmZWWxcuVKunTpcs5zFhUVkZ2dXWYTEXEWJpOJMb3CWXH/COZN6c/IHqFYrPDSD3u5+t9rGfbsT7z4/R6yCkrsXao0IsVmi71LqFSjDUIZGRmYzWZCQ0PL7A8NDSU1NRWA5ORkhg0bRnR0NMOGDeOuu+6iV69e5zznM888Q2BgoG2LjIys1/cgIuKIXF1MXNQ1lH/f2I9/jO/JsE4h+Hu5cSSzgFd+3M/4hDX8sj9Dg0+kSopLHTsIudm7gPoUGxtLYmJilY+fPXs2M2fOtD3Ozs5WGBIRp2UymbhpUBtuGtSGwhIz3+9M41/f7CYpI4/r31lHTGQQt1zQjou6tsTXs0n/OpFaKHHwFqFG+8kNCQnB1dW1XOfntLQ0wsLCanROT09PPD0966I8EZEmxcvdlSuiIxjSoTkvLdvLkk1HSDyUyV0fbaZVkDfPXNmLCzqGqA+RlLP9SJa9S6hUo7015uHhQb9+/Vi+fLltn8ViYfny5QwePLhW505ISKB79+4MGDCgtmWKiDQpIX6ePDWhFz8/EMftF7YnItCLI5kFTJq3njGvrGLuD3tZ+/sJe5cpDuT34445WuwMh24Rys3NZf/+/bbHSUlJJCYmEhwcTFRUFDNnzmTy5Mn079+f2NhY5s6dS15eHlOnTq3VdadPn8706dNtvc5FRKSslgFezB7TjTtHdOT57/awdPMRdqfmsDs1B9jHBR1DeGhUV3q11v9DnVlj6Efm0MPnV6xYwYgRI8rtnzx5MgsWLADgtdde47nnniM1NZWYmBheeeUVBg4cWCfX1/B5EZGqOZFbxFurDnDgeB4r9qRTYjZ+tfRsFcCoHmFcMyCSlv5a8d7Z7E/PJf7Fn22PHXH4vEMHIXtTEBIRqb5DJ/N5cdleliYe4cxvGHdXE3fGdeSO4R3w9nC1b4HSYC56YQUHzro1piDUSCQkJJCQkIDZbGbv3r0KQiIiNXAsq4DV+zL4aH0Km1IyAfB2d+Wiri35S2wkQzuoc3VT13bWV2UeKwg1MmoREhGpG59tPszz3+3lSGaBbV+XUH9uHdaO+G6hNPP1sGN1Ul/ODkIdWviy/K9xDXLd6vz+dujO0iIi0jRM6NOa8TGt2H4km083HebTjYfZk5bDA//dir+XG7dc0I4hHUIY0LYZJpNaiaThNNrh8yIi0riYTCZ6tQ7k8St6sPLBEUwb1o4OLXzJKSxl7g/7uObNtdz4n3XsTtXyRk3BidyiMo8d9faTgpCIiDS4Zr4e/O2y7nx/33Ceu7o3wzqF4OHqwpr9Jxj98irGvrqaJZsOY7Y46q9POZ/tRxtHoFUQqoAmVBQRaRiuLiYm9o/k/VsGsvyvw7msVzhWK2w7ksXMj7fQ6/Hv+Ntn2ziZV2zvUqWaCopLy+5w0EyrztKVUGdpEZGGl5ZdyKebDvPWygNk5v+x0r2nmwtX9m3NY2O74+WuIfiO7sN1yfzts+0E+biTmV9C+xBffrw/rkGuXZ3f32oREhERhxIa4MWdcR1Z9/DFzJvSn27hxi+yolILH61P4fJXV/PWyt85nlN0njOJPc39YR8AQd7ugMM2CGnUmIiIOCZPN1cu6hrKRV1DyS4sYWPyKe5bnMj+9Fye/no3//xmN3FdWvLEFT2IDPaxd7lyFqvVaguqB0/k27mayikIiYiIwwvwcmdEl5Z8f9+FfLMtlSWbj7DlUCY/7k7nx93pxEQGcWGnECb2j1QocgCfbjpi+z6+W0t+2JVux2oqp1tjFVBnaRERx9TS34vJQ9ry+fShfH33MAa1D8ZkgsRDmbzy437iX/yZf327m+1HsuxdqlO7/5Mttu9vuaA94LgLsKqzdCXUWVpExPEdySzgvV8OsmpfBjuP/TFku3+bZvSICGDq0Ha0DfG1Y4XO5+wZpT/9v8Fc9cZa2jb3YcUD5RdSrw+aWVpERJxGqyBvZo/pxiyrla+2HeOT3w6zZn8GvyWf4rfkU3y4LoW7LurEhD6tiGqu22b1LS270Pb9s1f1tn3vqK0uujUmIiJNgslk4vLeEbx7cywrHohjypC2AJRarLz0w15GvLCCuz7azNGz1juTuvfaj/tt31/ZtxXg2EumKAiJiEiT07qZD49f0YOtj1/KrNFd6Rzqh9li5X9bjjLs2Z+466PNrNp3XDNX14MzS6RM7NcaN9c/YoajdsTRrTEREWmyArzcuWN4B+4Y3oFNKad45utdbDh4iv9tOcr/thylY0s/xvQKZ2K/1hptVgdKzRb2p+cCcMOgNgA4+hq6ahGqgEaNiYg0PX2jmvHJHUP46u4LuH5gFAFebuxPz+WV5fuIe34Ft7//G2v2Z9i7zEbt+e/3ciq/hOa+HnQL9y/znNVBewlp1FglNGpMRKTpOpVXzPu/JrN6fwbrk07a9kdHBnHjwCjGxbTCw03tBdXR7x/LOJFXzJQhbXn8ih4AbE45xYTXfyEy2JtVD17UIHVo1JiIiMh5NPP14O6LO3H3xZ345fcMPt98lE83HWbLoUy2HMrkya92ce2ASG4a1Ea3zapg59FsTpxeHHf6iI52rqbqFIRERMTpDekQwpAOIdw5ogML16fw6cYjZOQW8dbKA7y96gAxkUH8Y1xPerYKtHepDmvemiQAYtsF08Lfs9zzjnr/SW1+IiIip7Vp7svs0d1Y9/DFvHJdH/q1aYbVCptTMrn81dWMmruSzzYf1mizP9l5NJv/bjwMwKB2wWWeMzl4b2kFIRERkT9xdTFxRXQEn/7fEFY9OILLeoVjMsHu1BzuW7yF/k8u49lvd5NfXGrvUh3CmFdW2b6fflHFt8UctUVIt8ZEREQqERnsQ8INfUnNKuTTTYd5Z9UBTuWX8PqK33lvbTKX9w5n8pC2dAt3zkE1Z4bLAzwwsguebq5lnnfs9iAFIRERkSoJC/Ri+oiO3H5he77fmcY/v9lNysl8Fm04xKINh+gS6s+0C9tzee9wvNxdz3/CJqDEbOG2938DwM3FxJ1xHexcUfXp1lgFNI+QiIici5urC2N6hbPi/jgW3zaI0T3DANiTlsP9n2zh4hd+5v21Bzl1egRVUzbr020cOJ4HwDuT+1fYH8jBuwhpHqHKaB4hERGpiv3pOXzwawpfbj1KRq4RgNxcTIzqGcYtF7QjJjLI4TsNV9eKPelMmb8BgEBvd7Y8dmmFx209nMkVr60hItCLX2Zf3CC1aR4hERGRBtSxpT+PX9GDh0Z1ZeH6FJZsOsyOo9l8ufUYX249RudQP64dEMWl3UObxJxE325P5f8+3AiAn6cbmx+95JzHmhy8l5CCkIiISB3x9nDllgvaccsF7dh5NJt3Vh3g6+3H2JuWyz++3Mk/vtzJsE4hzBjRkdh2wQ7fSlRcauHj3w6xOzWbQG93fDzceO67PbbnXV1MrHggDhcXx34flVEQEhERqQfdIwJ48doYHh/Xg88Tj7JwXQq7jmWzal8Gq/ZlEBbgxaU9QhnZI4zYdsG4uzpWt91th7MY+9rqcz7fq1Ug794cS7CvR5XO56j9cBSERERE6lGAlzs3DWrDTYPa8MvvGfz75wP8dvAkqdmFvLc2mffWJtPC35NB7ZszrFMIA9sF06a573nPeyqvmBN5RbQL8cO1Dltk8otLeeH7vfxnddI5j3l4TFemDWtfpRYtB2/0UhASERFpKGeW8igsMbNmfwbf70hj2a40jucU8b8tR/nflqMARAR6ERnsQ89WgYQHeuHt4YqXmyuZBSUkn8hjfdJJdqfmAEbQ+MuASG4a1JZu4f61ut2261g20xduso0Ea+nvybwpA+jZKpDswhJ83F1xdTHV6BqOOjRLQUhERKSBebm7cnG3UC7uFsoTJWbWJ520bb8ln+RoViFHswpZl3TyvOeyWuGj9Yf4aP0hooJ9uLx3OHFdWhIdGVhucsNz2Zh8iqe/3sWmlFO2wPLY2O5MGtzW1toU4OVe4/fryBSERERE7MjL3ZULO7fgws4tAMjML2ZvWi7rDpwg+WQ+e9Ny8PVww8UFmvt6EhrgSUxkM2LbBRPo7c6H65L5YVcavx08RcrJfF5f8Tuvr/gdgOjIIPpFNSM0wJNOoX6EBnhRVGohM7+YrIISNiafYtvhLLYczrLVc1nvcGaN6lrno9usDtpLSEFIRETEgQT5eBDbLpjYPy1eei5Th7Zj6tB25BWV8vW2Y6zcl8GPu9LIKzaz5VAmWw5lVuk87UJ8eeKKHrZAVlfUR6gRSkhIICEhAbPZbO9SREREqsTX042J/SOZ2D8Ss8XK0cwC1iedZNuRLHanZpOaVUhGbjGuLiaCfNwJ9feiZYAnfaOaMbpXGOGB3vVan6P2EdLM0pWoysyUVquV0tJShSZxCq6urri5uTn83Cci4jh2Hs1mzCuraOnvyfq/xTfINTWzdAMpLi7m2LFj5Ofn27sUkQbj4+NDeHg4Hh5VmztERMSRKQjVkMViISkpCVdXVyIiIvDw8NBfydKkWa1WiouLOX78OElJSXTq1AkXF8eaAE5EHJej3n5SEKqh4uJiLBYLkZGR+Pg0/nVjRKrC29sbd3d3kpOTKS4uxsvLy94liYiDc/Q2Av05V0v6i1icjT7zIlITjtojWf9HExERkXqjFiERERERB+0lpCDkhKZMmYLJZOKOO+4o99z06dMxmUxMmTKlSuc6ePAgJpOJxMTEctcYP3587YsVEZFGzYRjNwkpCDmpyMhIFi1aREFBgW1fYWEhCxcuJCoqyo6VlVVcXGzvEkREpAlTEHJSffv2JTIykiVLltj2LVmyhKioKPr06WPb9+2333LBBRcQFBRE8+bNufzyy/n9999tz7dr1w6APn36YDKZiIuL4/HHH+fdd9/l888/x2QyVilesWIFAIcOHeKaa64hKCiI4OBgxo0bx8GDB23nO9OS9NRTTxEREUGXLl3q9wchIiINwlE7S2v4fB2yWq0UlNhnhmlvd9dqz2N08803M3/+fG644QYA5s2bx9SpU22hBSAvL4+ZM2fSu3dvcnNzmTNnDhMmTCAxMREXFxfWr19PbGwsP/zwAz169MDDwwMPDw927dpFdnY28+fPByA4OJiSkhJGjhzJ4MGDWbVqFW5ubjz55JOMGjWKrVu32iboW758OQEBASxbtqxufjgiImI3jt5ZWkGoDhWUmOk+5zu7XHvn30fi41G9/5w33ngjs2fPJjk5GYA1a9awaNGiMkHoqquuKvOaefPm0aJFC3bu3EnPnj1p0cJYnK958+aEhYXZjvP29qaoqKjMvg8++ACLxcI777xjC23z588nKCiIFStWcOmllwLg6+vLO++8o5mLRUSaEAdtEFIQcmYtWrTgsssuY8GCBVitVi677DJCQkLKHLNv3z7mzJnDunXryMjIwGKxAJCSkkLPnj2rdb0tW7awf/9+/P39y+wvLCwsc7utV69eCkEiIk2EgzcIKQhVpKarz3u7u7Lz7yPrqarzX7smbr75ZmbMmAEY7/vPxo4dS5s2bXj77beJiIjAYrHQs2fPGnVizs3NpV+/fnz44YflnjvTsgRGi5CIiDQtjrrGu4JQBaZPn8706dNtq9dWlclkqvbtKXsbNWoUxcXFmEwmRo4sG+JOnDjBnj17ePvttxk2bBgAq1evLnPMmZabP4dGDw+Pcvv69u3L4sWLadmy5XlXAxYRkabB0fsIadSYk3N1dWXXrl3s3LkTV9eyrUrNmjWjefPmvPXWW+zfv58ff/yRmTNnljmmZcuWeHt78+2335KWlkZWVhYAbdu2ZevWrezZs4eMjAxKSkq44YYbCAkJYdy4caxatYqkpCRWrFjB3XffzeHDhxvsPYuISMNzzPYgBSEBAgICKmyhcXFxYdGiRWzcuJGePXty33338dxzz5U5xs3NjVdeeYU333yTiIgIxo0bB8C0adPo0qUL/fv3p0WLFqxZswYfHx9WrlxJVFQUV155Jd26deOWW26hsLBQLUQiIk2WYzcJmayOetPOAZy5NZaVlVXuF3VhYSFJSUm0a9dOK3CLU9FnX0SqY396LvEv/kyQjzuJcy5tkGtW9vv7z9QiJCIiIvXOUZtdFIRERESk3qiztIiIiDg9R+2JoyAkIiIi9cbBG4QUhERERKT+OWZ7kIKQiIiI1KPqLgje0BSEREREpP45aJOQgpCIiIjUG8duD1IQEhERESemICQ1FhcXx7333mvvMkREpBFw0DtjCkLOaMqUKYwfP77c/hUrVmAymcjMzGzwmkREpGly8L7SCkIiIiJS/zShojQqJ06c4LrrrqNVq1b4+PjQq1cvPvroo0pfc+rUKSZNmkSzZs3w8fFh9OjR7Nu3z/Z8cnIyY8eOpVmzZvj6+tKjRw++/vrr+n4rIiJiRyYH7y7tZu8C6tOhQ4e46aabSE9Px83NjUcffZSJEyfW3wWtVijJr7/zV8bdp07bHwsLC+nXrx8PPfQQAQEBfPXVV9x000106NCB2NjYCl8zZcoU9u3bxxdffEFAQAAPPfQQY8aMYefOnbi7uzN9+nSKi4tZuXIlvr6+7Ny5Ez8/vzqrWUREHJdjtgc18SDk5ubG3LlziYmJITU1lX79+jFmzBh8fX3r54Il+fB0RP2c+3wePgoeVX9fX375ZbkQYjabbd+3atWK+++/3/b4rrvu4rvvvuPjjz+uMAidCUBr1qxhyJAhAHz44YdERkaydOlSJk6cSEpKCldddRW9evUCoH379tV6iyIi0vg4eh+hJh2EwsPDCQ8PByAsLIyQkBBOnjxZf0GoERkxYgRvvPFGmX3r1q3jxhtvBIxQ9PTTT/Pxxx9z5MgRiouLKSoqwsfHp8Lz7dq1Czc3NwYOHGjb17x5c7p06cKuXbsAuPvuu/m///s/vv/+e+Lj47nqqqvo3bt3Pb1DERFxJA7aRcixg9DKlSt57rnn2LhxI8eOHeOzzz4rN9opISGB5557jtTUVKKjo3n11VcrbLHYuHEjZrOZyMjI+ivY3cdombEH94oDyrn4+vrSsWPHMvsOHz5s+/65557j5ZdfZu7cufTq1QtfX1/uvfdeiouLa1zirbfeysiRI/nqq6/4/vvveeaZZ3jhhRe46667anxOERGR2nDoztJ5eXlER0eTkJBQ4fOLFy9m5syZPPbYY2zatIno6GhGjhxJenp6meNOnjzJpEmTeOutt+q3YJPJuD1lj62O2x7XrFnDuHHjuPHGG4mOjqZ9+/bs3bv3nMd369aN0tJS1q1bZ9t34sQJ9uzZQ/fu3W37IiMjueOOO1iyZAl//etfefvtt+u0bhERkepw6CA0evRonnzySSZMmFDh8y+++CLTpk1j6tSpdO/enX//+9/4+Pgwb9482zFFRUWMHz+eWbNm2fqunEtRURHZ2dllNmfVqVMnli1bxi+//MKuXbu4/fbbSUtLq/T4cePGMW3aNFavXs2WLVu48cYbadWqFePGjQPg3nvv5bvvviMpKYlNmzbx008/0a1bt4Z6SyIiYkdWB+0u7dBBqDLFxcVs3LiR+Ph42z4XFxfi4+NZu3YtYMxZMGXKFC666CJuuumm857zmWeeITAw0LbV6200B/fII4/Qt29fRo4cSVxcHGFhYRVOwni2+fPn069fPy6//HIGDx6M1Wrl66+/xt3dHTD6HU2fPp1u3boxatQoOnfuzOuvv94A70ZEROxFnaXrSUZGBmazmdDQ0DL7Q0ND2b17N2Dc3lm8eDG9e/dm6dKlALz//vu2UUt/Nnv2bGbOnGl7nJ2d3STD0IIFCyrcHxcXV2bCqzM/s3NZsWJFmcfNmjXjvffeO+fxr776alVLFBGRJkadpe3gggsuwGKxVPl4T09PPD0967EiERER52Jy8CahRntrLCQkBFdX13L9VtLS0ggLC6vVuRMSEujevTsDBgyo1XlERETE4KANQo03CHl4eNCvXz+WL19u22exWFi+fDmDBw+u1bmnT5/Ozp072bBhQ23LFBERcWqO3R7k4LfGcnNz2b9/v+1xUlISiYmJBAcHExUVxcyZM5k8eTL9+/cnNjaWuXPnkpeXx9SpU+1YtYiIiJTjoE1CDh2EfvvtN0aMGGF7fKYj8+TJk1mwYAHXXnstx48fZ86cOaSmphITE8O3335brgO1iIiI2IeDdxFy7CD051FMFZkxYwYzZsyo0+smJCSQkJBQZu0tERERaXoabR+h+qQ+QiIiInVLEyqKiIiI0zE5eHdpBSERERGpd446oaKCkDiFuLg47r333hq9dsGCBQQFBdVpPfXt8ccfJyYmxt5liIg4fGdpBaEKNPUJFadMmYLJZCq3jRo1yt6llVGb8OJoDh48iMlkIjExsUGud//995eZY0tExN4ctEHIsUeN2cv06dOZPn062dnZBAYG2rucejFq1Cjmz59fZp+WF2k6/Pz88PPzs3cZIiIO3kNILUJOy9PTk7CwsDJbs2bNAGMxVQ8PD1atWmU7/tlnn6Vly5a2JU3i4uJsUxcEBgYSEhLCo48+Wma6g6KiIu6//35atWqFr68vAwcOLLdQ65o1a4iLi8PHx4dmzZoxcuRITp06xZQpU/j55595+eWXbS1WBw8eBGD79u2MHj0aPz8/QkNDuemmm8jIyLCdMy8vj0mTJuHn50d4eDgvvPDCeX8eW7ZsYcSIEfj7+xMQEEC/fv347bffyhzz3Xff0a1bN/z8/Bg1ahTHjh2zPWexWPj73/9O69at8fT0tM1pdUa7du0A6NOnDyaTibi4uHPWEhcXx913382DDz5IcHAwYWFhPP7442WOSUlJYdy4cfj5+REQEMA111xTZrmZP98aW7FiBbGxsfj6+hIUFMTQoUNJTk62Pf/555/Tt29fvLy8aN++PU888QSlpaXn/bmJiFTV+abDsRcFoTpktVrJL8m3y1aXH7Azt6RuuukmsrKy2Lx5M48++ijvvPNOmckq3333Xdzc3Fi/fj0vv/wyL774Iu+8847t+RkzZrB27VoWLVrE1q1bmThxIqNGjWLfvn0AJCYmcvHFF9O9e3fWrl3L6tWrGTt2LGazmZdffpnBgwczbdo0jh07xrFjx4iMjCQzM5OLLrqIPn368Ntvv/Htt9+SlpbGNddcY7vuAw88wM8//8znn3/O999/z4oVK9i0aVOl7/mGG26gdevWbNiwgY0bNzJr1izc3d1tz+fn5/P888/z/vvvs3LlSlJSUrj//vttz7/88su88MILPP/882zdupWRI0dyxRVX2N7r+vXrAfjhhx84duwYS5YsqbSed999F19fX9atW8ezzz7L3//+d5YtWwYYoWvcuHGcPHmSn3/+mWXLlnHgwAGuvfbaCs9VWlrK+PHjGT58OFu3bmXt2rXcdttttoUQV61axaRJk7jnnnvYuXMnb775JgsWLOCpp56qtEYRkSpx8CYhk9VRI5oDOHNrLCsri4CAgDLPFRYWkpSURLt27fDy8gIgvySfgQsH2qNU1l2/Dh93nyodO2XKFD744ANb3Wc8/PDDPPzwwwAUFxczcOBAOnfuzPbt2xk6dChvvfWW7di4uDjS09PZsWOH7RfqrFmz+OKLL9i5cycpKSm0b9+elJQUIiIibK+Lj48nNjaWp59+muuvv56UlBRWr15dYZ1xcXHExMQwd+5c274nn3ySVatW8d1339n2HT58mMjISPbs2UNERATNmzfngw8+YOLEiQCcPHmS1q1bc9ttt5U519kCAgJ49dVXmTx5crnnFixYwNSpU9m/fz8dOnQA4PXXX+fvf/87qampALRq1Yrp06fbfn4AsbGxDBgwgISEBA4ePEi7du3YvHnzeTsxx8XFYTaby7TIxcbGctFFF/HPf/6TZcuWMXr0aJKSkoiMjARg586d9OjRg/Xr1zNgwAAef/xxli5dSmJiIidPnqR58+asWLGC4cOHl7tefHw8F198MbNnz7bt++CDD3jwwQc5evRoueMr+uyLiJxLek4hsU8tx8UEB565rEGuWdnv7z9TH6EKOMPM0iNGjOCNN94osy84ONj2vYeHBx9++CG9e/emTZs2vPTSS+XOMWjQIFsIAhg8eDAvvPACZrOZbdu2YTab6dy5c5nXFBUV0bx5c8BoEToTVqpqy5Yt/PTTTxX2f/n9998pKCiwhbiz31eXLl0qPe/MmTO59dZbef/994mPj2fixIm20APg4+NT5nF4eDjp6emA8Q/u6NGjDB06tMw5hw4dypYtW855zVWrVjF69Gjb4zfffJMbbrgBgN69e5c59uzr7dq1i8jISFsIAujevTtBQUHs2rWrXCf/4OBgpkyZwsiRI7nkkkuIj4/nmmuuITw8HDB+pmvWrCnTAmQ2myksLCQ/Px8fn6oFbBGRyjhqq4uCUAVq2lna282bddevq8fKKr92dfj6+tKxY8dKj/nll18Ao0Xl5MmT+Pr6Vvn8ubm5uLq6snHjRlxdXcs8dybEeHtXr+Yz5x07diz/+te/yj0XHh5eZpHe6nj88ce5/vrr+eqrr/jmm2947LHHWLRoERMmTAAoc5sMwGQy1fp2ZP/+/cuMIjv7tmNF17NYLDW+1vz587n77rv59ttvWbx4MY888gjLli1j0KBB5Obm8sQTT3DllVeWe51afESkthx9QkUFoTpkMpmqfHvK0f3+++/cd999vP322yxevJjJkyfzww8/4OLyR7eydevKhr5ff/2VTp064erqSp8+fTCbzaSnpzNs2LAKr9G7d2+WL1/OE088UeHzHh4e5Vrl+vbty6effkrbtm1xcyv/8e3QoQPu7u6sW7eOqKgoAE6dOsXevXsrvC10ts6dO9O5c2fuu+8+rrvuOubPn28LQpUJCAggIiKCNWvWlLnGmjVriI2Ntb0XoMz78fb2Pm8YrUi3bt04dOgQhw4dKnNrLDMzk+7du5/zdX369KFPnz7Mnj2bwYMHs3DhQgYNGkTfvn3Zs2dPjWoREakqR+2Io87STqqoqIjU1NQy25mRV2azmRtvvJGRI0cydepU5s+fz9atW8uNvkpJSWHmzJns2bOHjz76iFdffZV77rkHMELFDTfcwKRJk1iyZAlJSUmsX7+eZ555hq+++gqA2bNns2HDBu688062bt3K7t27eeONN2x1tG3blnXr1nHw4EEyMjKwWCxMnz6dkydPct1117FhwwZ+//13vvvuO6ZOnYrZbMbPz49bbrmFBx54gB9//JHt27czZcqUMgHuzwoKCpgxYwYrVqwgOTmZNWvWsGHDBrp161bln+cDDzzAv/71LxYvXsyePXuYNWsWiYmJtp9Hy5Yt8fb2tnXuzsrKqvp/rD+Jj4+nV69e3HDDDWzatIn169czadIkhg8fTv/+/csdn5SUxOzZs1m7di3Jycl8//337Nu3z/b+5syZw3vvvccTTzzBjh072LVrF4sWLeKRRx6pcY0iImc4+oSKWOWcsrKyrIA1Kyur3HMFBQXWnTt3WgsKCuxQWe1MnjzZinG7tszWpUsXq9VqtT7xxBPW8PBwa0ZGhu01n376qdXDw8OamJhotVqt1uHDh1vvvPNO6x133GENCAiwNmvWzPrwww9bLRaL7TXFxcXWOXPmWNu2bWt1d3e3hoeHWydMmGDdunWr7ZgVK1ZYhwwZYvX09LQGBQVZR44caT116pTVarVa9+zZYx00aJDV29vbCliTkpKsVqvVunfvXuuECROsQUFBVm9vb2vXrl2t9957r+3aOTk51htvvNHq4+NjDQ0NtT777LPW4cOHW++5554Kfx5FRUXWv/zlL9bIyEirh4eHNSIiwjpjxgzbf9v58+dbAwMDy7zms88+s579z8dsNlsff/xxa6tWrazu7u7W6Oho6zfffFPmNW+//bY1MjLS6uLiYh0+fPg5//tUVOu4ceOskydPtj1OTk62XnHFFVZfX1+rv7+/deLEidbU1FTb84899pg1OjraarVarampqdbx48dbw8PDrR4eHtY2bdpY58yZYzWbzbbjv/32W+uQIUOs3t7e1oCAAGtsbKz1rbfeqrC+xvzZF5GGl5lXbL36jTXWq99Y02DXrOz3959p1FglqjtqzJlUNKJLnIOzf/ZFxPFVZ9SYbo2JiIiI01IQqkBTX2tMREREDLo1VgndGhMpT599EXF0ujUmIiIiUgUKQiIiIuK0FIRqSXcWxdnoMy8iTYmCUA2dWQIhPz/fzpWINKwzn/k/LwMiItIYaYmNGnJ1dSUoKMi2EKaPj0+ZBUhFmhqr1Up+fj7p6ekEBQWVW0NORKQxUhCqQFVXnw8LCwOwhSERZxAUFGT77IuINHYaPl+Jqg6/M5vNlJSUNGBlIvbh7u6uliARcXjVGT6vFqE64Orqql8OIiIijZA6S4uIiIjTUhASERERp6UgJCIiIk5LfYQqcaYfeXZ2tp0rERERkao683u7KuPBFIQqkZOTA0BkZKSdKxEREZHqysnJITAwsNJjNHy+EhaLhaNHj+Lv72+bLHHAgAFs2LChSq+vyrGVHZOdnU1kZCSHDh067/C/xqY6P8fGcu26OG9tzlHd11b1eH2Oz60pfo7r6tw1PYejfo6h6X6W7fk5rq/rW61W+vXrx969e3FxqbwXkFqEKuHi4kLr1q3L7HN1da3yP4CqHFuVYwICAprUPzqo3s+xsVy7Ls5bm3NU97VVPV6f43Nrip/jujp3Tc/h6J9jaHqfZXt+juvz+h4eHucNQaDO0tU2ffr0Oj22OudrSuz5vuvr2nVx3tqco7qvrerx+hyfW1P8HNfVuWt6Dn2OG56937e9/5+sW2MOrDozY4o4Kn2OpanQZ7lpUouQA/P09OSxxx7D09PT3qWI1Jg+x9JU6LPcNKlFSERERJyWWoRERETEaSkIiYiIiNNSEBIRERGnpSAkIiIiTktBqAk4dOgQcXFxdO/end69e/PJJ5/YuySRGpswYQLNmjXj6quvtncpIlX25Zdf0qVLFzp16sQ777xj73KkGjRqrAk4duwYaWlpxMTEkJqaaptW3NfX196liVTbihUryMnJ4d133+W///2vvcsROa/S0lK6d+/OTz/9RGBgIP369eOXX36hefPm9i5NqkAtQk1AeHg4MTExAISFhRESEsLJkyftW5RIDcXFxeHv72/vMkSqbP369fTo0YNWrVrh5+fH6NGj+f777+1dllSRglADWLlyJWPHjiUiIgKTycTSpUvLHZOQkEDbtm3x8vJi4MCBrF+/vkbX2rhxI2azmcjIyFpWLVJeQ36WRRpKbT/XR48epVWrVrbHrVq14siRIw1RutQBBaEGkJeXR3R0NAkJCRU+v3jxYmbOnMljjz3Gpk2biI6OZuTIkaSnp9uOiYmJoWfPnuW2o0eP2o45efIkkyZN4q233qr39yTOqaE+yyINqS4+19KIWaVBAdbPPvuszL7Y2Fjr9OnTbY/NZrM1IiLC+swzz1T5vIWFhdZhw4ZZ33vvvboqVaRS9fVZtlqt1p9++sl61VVX1UWZItVSk8/1mjVrrOPHj7c9f88991g//PDDBqlXak8tQnZWXFzMxo0biY+Pt+1zcXEhPj6etWvXVukcVquVKVOmcNFFF3HTTTfVV6kilaqLz7KIo6nK5zo2Npbt27dz5MgRcnNz+eabbxg5cqS9SpZqUhCys4yMDMxmM6GhoWX2h4aGkpqaWqVzrFmzhsWLF7N06VJiYmKIiYlh27Zt9VGuyDnVxWcZID4+nokTJ/L111/TunVrhSixq6p8rt3c3HjhhRcYMWIEMTEx/PWvf9WIsUbEzd4FSO1dcMEFWCwWe5chUid++OEHe5cgUm1XXHEFV1xxhb3LkBpQi5CdhYSE4OrqSlpaWpn9aWlphIWF2akqkerTZ1maIn2umz4FITvz8PCgX79+LF++3LbPYrGwfPlyBg8ebMfKRKpHn2VpivS5bvp0a6wB5Obmsn//ftvjpKQkEhMTCQ4OJioqipkzZzJ58mT69+9PbGwsc+fOJS8vj6lTp9qxapHy9FmWpkifaydn72FrzuCnn36yAuW2yZMn24559dVXrVFRUVYPDw9rbGys9ddff7VfwSLnoM+yNEX6XDs3rTUmIiIiTkt9hERERMRpKQiJiIiI01IQEhEREaelICQiIiJOS0FIREREnJaCkIiIiDgtBSERERFxWgpCIiIi4rQUhERERMRpKQiJiIiI01IQEhFpABMmTKBZs2ZcffXV9i5FRM6iICQi0gDuuece3nvvPXuXISJ/oiAkIg3q/vvvZ/z48VU6Ni4uDpPJhMlkIjExsUbncBRxcXH4+/tX+NyUKVNs73Pp0qUNW5iIk1MQEpEGlZiYSExMTJWPnzZtGseOHaNnz55lztG7d+9zvuZMsLjjjjvKPTd9+nRMJhNTpkypTtn16uWXX+bYsWP2LkPEKSkIiUiD2rJlS7WCkI+PD2FhYbi5uZU5R3R0dKWvi4yMZNGiRRQUFNj2FRYWsnDhQqKioqpd9/nExMTQs2fPctvRo0fP+9rAwEDCwsLqvCYROT8FIRFpMIcPHyYjI8MWhDIzMxk7diwXXHABqamp1ToHwCWXXIKPjw9dunRh3bp1ZY7r27cvkZGRLFmyxLZvyZIlREVF0adPnzLHxsXFMWPGDGbMmEFgYCAhISE8+uijWK1W2zEWi4Vnn32Wjh074unpSVRUFE899ZTt+cTERLZv315ui4iIqNbPSEQaloKQiDSYxMREgoKCaNu2Ldu2bWPAgAG0atWKn376qcotImf6CiUkJPDwww+zZcsWoqKimDVrVrljb775ZubPn297PG/ePKZOnVrhed99913c3NxYv349L7/8Mi+++CLvvPOO7fnZs2fzz3/+k0cffZSdO3eycOFCQkNDq/HuRcQRuZ3/EBGRupGYmEh0dDQLFy5kxowZ/Otf/2LatGnVPkdwcDAff/wxISEhAFxxxRW8+eab5Y698cYbmT17NsnJyQCsWbOGRYsWsWLFinLHRkZG8tJLL2EymejSpQvbtm3jpZdeYtq0aeTk5PDyyy/z2muvMXnyZAA6dOjABRdcUOW64+Pj2bJlC3l5ebRu3ZpPPvmEwYMHV+u9i0jdUxASkQaTmJjI1q1bmTFjBl999VWNgkBiYiLjxo2zhSCApKQkOnbsWO7YFi1acNlll7FgwQKsViuXXXZZmdedbdCgQZhMJtvjwYMH88ILL2A2m9m1axdFRUVcfPHF1a73jB9++KHGrxWR+qNbYyLSYBITE7nyyispLCwkMzOzxucYNGhQuX3n6oB98803s2DBAt59911uvvnmGl3T29u7Rq8TEcenICQiDSInJ4cDBw4wffp0XnvtNf7yl7+wY8eOGp3jz52dKwtCo0aNori4mJKSEkaOHHnOc/+5s/Wvv/5Kp06dcHV1pVOnTnh7e7N8+fJq1Ssijk+3xkSkQWzZsgVXV1e6d+9Onz592L59O2PHjmX9+vXnvF11rnP06tXLti85OZlTp06dMwi5urqya9cu2/fnkpKSwsyZM7n99tvZtGkTr776Ki+88AIAXl5ePPTQQzz44IN4eHgwdOhQjh8/zo4dO7jllluq+BMQEUekICQiDSIxMZGuXbvi6ekJwHPPPceuXbu48sor+eGHH/Dw8KjSObp06YKXl5dt3+bNm20j0c4lICDgvOeeNGkSBQUFxMbG4urqyj333MNtt91me/7RRx/Fzc2NOXPmcPToUcLDwyucsFFEGheT9eyJMkREHEhcXBwxMTHMnTu3SVznfEwmE5999lmjWz5EpDFTHyERcWivv/46fn5+bNu2zd6l1Js77rgDPz8/e5ch4pTUIiQiDuvIkSO2JTKioqKqdPusJuzdIpSenk52djYA4eHh+Pr62qUOEWekICQiIiJOS7fGRERExGkpCImIiIjTUhASERERp6UgJCIiIk5LQUhEREScloKQiIiIOC0FIREREXFaCkIiIiLitBSERERExGkpCImIiIjTUhASERERp6UgJCIiIk5LQUhERESc1v8DRWPHxTrxee4AAAAASUVORK5CYII=\n"},"metadata":{}}]},{"metadata":{},"cell_type":"markdown","source":"As can be seen, halos are more strongly clustered than matter, as expected as we are taking galaxy clusters with masses above 1e14 Msun/h. On small scales, the power spectrum of halos saturates at the expected shot-noise level. On the smallest scales we have, the power spectrum of both halos and matter is affected by aliasing so should not be trusted."},{"metadata":{},"cell_type":"markdown","source":"### Finally, lets show an example of how to compute marked power spectra"},{"metadata":{},"cell_type":"markdown","source":"Our goal is to compute the power spectrum of halos but instead of giving each halo the same weight, we want to weight each halo by the overdensity of neutrinos from the cosmic neutrino background."},{"metadata":{"trusted":true},"cell_type":"code","source":"snapshot = '/home/jovyan/Data/Snapshots/Mnu_ppp/261/snapdir_004/snap_004' #location of the snapshot\nsnapdir  = '/home/jovyan/Data/Halos/FoF/Mnu_ppp/261/'                     #folder hosting the catalogue\nsnapnum  = 4                                                              #number of the catalog (4-->z=0, 3-->z=0.5, 2-->z=1, 1-->z=2, 0-->z=3)\n\n# read header\nheader   = readgadget.header(snapshot)\nBoxSize  = header.boxsize/1e3  #Mpc/h\nNall     = header.nall         #Total number of particles\nMasses   = header.massarr*1e10 #Masses of the particles in Msun/h\nOmega_m  = header.omega_m      #value of Omega_m\nOmega_l  = header.omega_l      #value of Omega_l\nh        = header.hubble       #value of h\nredshift = header.redshift     #redshift of the snapshot\nHubble   = 100.0*np.sqrt(Omega_m*(1.0+redshift)**3+Omega_l)#Value of H(z) in km/s/(Mpc/h)\n\n# read the neutrino positions\npos_n = readgadget.read_block(snapshot, \"POS \", [2])/1e3 #positions in Mpc/h","execution_count":14,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now, lets compute the density field of neutrinos"},{"metadata":{"trusted":true},"cell_type":"code","source":"grid = 512\n\n# define the matrix that will contain the neutrino density field\ndelta_n = np.zeros((grid,grid,grid), dtype=np.float32)\n\n# compute the neutrino density field\nMASL.MA(pos_n, delta_n, BoxSize, MAS, verbose=verbose)\n\n# compute the overdensity field\ndelta_n /= np.mean(delta_n, dtype=np.float64);  delta_n -= 1.0\n\n# print some information\nprint('%.3f < delta < %.3f'%(np.min(delta_n), np.max(delta_n)))\nprint('< delta > = %.3f'%np.mean(delta_n))","execution_count":15,"outputs":[{"output_type":"stream","text":"\nUsing CIC mass assignment scheme\nTime taken = 23.728 seconds\n\n-1.000 < delta < 51.771\n< delta > = 0.000\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"Now, lets read the halo catalog"},{"metadata":{"trusted":true},"cell_type":"code","source":"# read the halo catalogue\nFoF = readfof.FoF_catalog(snapdir, snapnum, long_ids=False,\n                          swap=False, SFR=False, read_IDs=False)\n\n# get the properties of the halos\npos_h  = FoF.GroupPos/1e3            #Halo positions in Mpc/h\nvel_h  = FoF.GroupVel*(1.0+redshift) #Halo peculiar velocities in km/s\nmass_h = FoF.GroupMass*1e10          #Halo masses in Msun/h\nNp_h   = FoF.GroupLen                #Number of CDM particles in the halo. Even in simulations with massive neutrinos, this will be just the number of CDM particles","execution_count":16,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We need to compute the overdensity of neutrinos in the location of the dark matter halos. For this we do"},{"metadata":{"trusted":true},"cell_type":"code","source":"# definte the array hosting the neutrino overdensities\ndelta_n_h = np.zeros(pos_h.shape[0], dtype=np.float32)\n\n# interpolate to find the neutrino overdensity in the positions of the halos\nMASL.CIC_interp(delta_n, BoxSize, pos_h, delta_n_h)","execution_count":17,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now, we can construct a density field weigthing each halo by its neutrino overdensity"},{"metadata":{"trusted":true},"cell_type":"code","source":"# matrix that will host the density field\ndelta = np.zeros((grid,grid,grid), dtype=np.float32)\n\n# compute the \"marked\" halo density field\nMASL.MA(pos_h, delta, BoxSize, MAS, W=delta_n_h, verbose=verbose)\n\n# print some information about the density field\nprint('%.3f < delta < %.3f'%(np.min(delta), np.max(delta)))\nprint('< delta > = %.3f'%np.mean(delta))","execution_count":18,"outputs":[{"output_type":"stream","text":"\nUsing CIC mass assignment scheme with weights\nTime taken = 0.520 seconds\n\n-0.808 < delta < 21.207\n< delta > = 0.002\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"We now compute the power spectrum of this field. This is equivalent to say that we are computing the marked power spectrum of the halos where the mark is the neutrino overdensity"},{"metadata":{"trusted":true},"cell_type":"code","source":"# compute power spectrum\naxis = 0 #we are working in real-space, so this value doesnt matter\nPk_h = PKL.Pk(delta, BoxSize, axis, MAS, threads, verbose)\n\n# Pk is a python class containing the 1D, 2D and 3D power spectra, that can be retrieved as\n\n# 1D P(k)\nk1D_h      = Pk_h.k1D\nPk1D_h     = Pk_h.Pk1D\nNmodes1D_h = Pk_h.Nmodes1D\n\n# 2D P(k)\nkpar_h     = Pk_h.kpar\nkper_h     = Pk_h.kper\nPk2D_h     = Pk_h.Pk2D\nNmodes2D_h = Pk_h.Nmodes2D\n\n# 3D P(k)\nk_h       = Pk_h.k3D\nPk0_h     = Pk_h.Pk[:,0] #monopole\nPk2_h     = Pk_h.Pk[:,1] #quadrupole\nPk4_h     = Pk_h.Pk[:,2] #hexadecapole\nPkphase_h = Pk_h.Pkphase #power spectrum of the phases\nNmodes_h  = Pk_h.Nmodes3D","execution_count":19,"outputs":[{"output_type":"stream","text":"\nComputing power spectrum of the field...\nTime to complete loop = 7.23\nTime taken = 12.47 seconds\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"Now lets see how this marked power spectrum looks like:"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.xscale('log')\nplt.yscale('log')\nplt.xlabel(r'$k~[h{\\rm Mpc}^{-1}]$')\nplt.ylabel(r'$P(k)~[(h^{-1}{\\rm Mpc})^3]$')\nplt.plot(k, Pk0_h)\nplt.show()","execution_count":20,"outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 640x480 with 1 Axes>","image/png":"iVBORw0KGgoAAAANSUhEUgAAAkoAAAG9CAYAAADqXFmlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABRhUlEQVR4nO3dd3RUZeLG8e9Meg8kJJACoYQSSgIBJIiKK4gNFAQbKoIiKpYV17au4vpb17bi7moEBSk2UBRZewEpSg0loRdpCQlJgEAapM78/ghGEQLJZJI7M3k+5+QcMnNz5wlnIE/e+973NVmtVisiIiIicgaz0QFEREREHJWKkoiIiEgNVJREREREaqCiJCIiIlIDFSURERGRGqgoiYiIiNRARUlERESkBipKIiIiIjVwNzqAM7NYLGRlZREQEIDJZDI6joiIiNSC1WqlsLCQiIgIzOZzjxmpKNVDVlYW0dHRRscQERERG2RkZBAVFXXOY1SU6iEgIACo+osODAw0OI2IiIjURkFBAdHR0dU/x89FRakefr3cFhgYqKIkIiLiZGozbUaTuUVERERqoKIkIiIiUgMVJREREZEaqCiJiIiI1EBFSURERKQGKkoiIiIiNVBREhEREamBipKIiIhIDVSURERERGqgoiQiIiJSAxUlERERkRqoKImIiIjUQEVJHJ7VauVQ/kl+2n2YE2UVRscREZEmxN3oACJ/dLSolE0H8099HCftYD5HikoB6Nk6mLnj++Ht4WZwShERaQpUlMRQBSXlbDmYT9rBfDZnHictI5/M4yfPOM7NbMLNbGJj+nEe+TiN12/uidlsMiCxiIg0JSpK0mhOllWy7VA+aRlVI0WbMvPZe7j4rMe2a+FHfFQwPaKC6BEVRFyrINIOHue2d9bw1eZDtA7x5fErOjfydyAiIk2NipI0iLIKC7tyCkk7eJxNGfmkHTzO7twiKi3WM46NDPYhPjqIHqeKUbfIIAK9Pc44rl+7EF4c0YNH5qcxdeke2jT35aa+rRvj2xERkSZKRUnqrdJiZc/hotPmFG0/VEBZheWMY0P9vYiPOlWKooPoERlEiL9XrV/r+sQoDuSd4L+Ld/O3hVuIaubLgNhQe347IiIi1VSUxCYWi5WZK/bx/bYctmbmU1xWecYxgd7u1aNEPaKCiY8OomWgNyZT/eYWPTwolgNHi/lfahb3vr+eT+/rT8fwgHqdU0RE5GxUlKTOSsorefijVL7Zkl39mI+HG90jg+h+ak5RfFQwbUJ8612KzsZkMvHyyB5kHT9Jyv5jjJ2VwsKJF9IioPYjUyIiIrVhslqtZ04akVopKCggKCiI/Px8AgMDjY7TKI4WlXLXu+vYmH4cTzczjw7pxMUdW9AhzB+3Rr4L7VhxGSOmrmTfkWLio4OZN74fPp5aNkBERM6tLj+/teCk1Nqew0UMf3MlG9OPE+Tjwbt39mX8xe3o1DKg0UsSQDM/T2be0YdgXw/SMo7z8EepWM4yWVxERMRWKkpSK2v35THizZWk550gurkPn97bn37tQoyORdtQP96+rTeebma+3ZrNS9/uMDqSiIi4EBUlOa//pWZy64w15J8sJz46mM/uu5AOYf5Gx6rWt21zXhnVA4C3lu/lwzXpBicSERFXoaIkNbJarSQv+YWH5qVSVmlhSNdw5o3vR2gdbudvLNcmRPLwoI4APP2/LSzfddjgRCIi4gpUlOSsyistPPHpZl75bicAdw1oy5ujEx16svSDl3VgRM9IKi1W7vtgAzuzC42OJCIiTk5FSc5QWFLOuNkpfLQuA7MJ/j6sK3+7Js6QCdt1YTKZeOH67lzQtjlFpRWMm51CbkGJ0bFERMSJqSjJabKOn2TUtFX8tPsIPh5uvH1bb8b0jzE6Vq15ubvx1m2JtAv1I/P4Se56dx0nyiqMjiUiIk5KRUmqbc3KZ/ibK9iRXUiLAC8+npDEoLhwo2PVWbCvJ7PG9qG5nyebDubz0LzUs+4xJyIicj4qSgLAkh253DBtFTkFpXQM9+ez+/rTPSrI6Fg2axPix9u3JeLpbuaHbTm88PV2oyOJiIgTUlES3l99gDvnpFBcVkn/9iHMv6c/Uc18jY5Vb71jmvOvUfEAzPh5H++tPmBwIhERcTYqSk2YxWLlha+387eFW7BYYWRiFLPH9iXIx8PoaHYzLD6CR4d0AmDy/7awZGeuwYlERMSZNPmi9OWXX9KpUydiY2OZMWOG0XEaTUl5JQ/M3chby/cCMGlwR14Z2QNPd9d7S9w3sD0jE6OwWOH+DzawLavA6EgiIuIkmvSmuBUVFcTFxbFkyRKCgoJITExk5cqVhITUbmsOZ90U92hRKePfXceG9ON4uJl4eWQPhveMMjpWgyqrsDBm5lpW7T1KqyBvFk68kPBAb6NjiYiIAbQpbi2tXbuWrl27EhkZib+/P1deeSXff/+90bEa1N7DRYyYupIN6ccJ9Hbn3XEXuHxJAvB0NzPt1kTat/DjUH4J42anUFyqZQNEROTcnLooLV++nKFDhxIREYHJZGLhwoVnHJOcnExMTAze3t5ccMEFrF27tvq5rKwsIiMjqz+PjIwkMzOzMaIbImV/HiOmruTA0RNENfNhwX39SWpv/Ma2jSXI14NZd/QlxM+TrVkFPDRvo5YNEBGRc3LqolRcXEx8fDzJyclnff6jjz5i0qRJTJ48mQ0bNhAfH8+QIUPIzW16E3q/SMti9PQ1HD9RTnxU0KmNbQOMjtXoWof4Mn1Mb7zczSzanss/vtpmdCQREXFgTl2UrrzySv7xj38wfPjwsz4/ZcoUxo8fz9ixY4mLi2PatGn4+voyc+ZMACIiIk4bQcrMzCQiIqLG1ystLaWgoOC0D0dntVp5c+kvPDB3I2WVFi6PC2fe3Um0CHC8jW0bS6/WzZhyQwIAs1bs55P1B40NJCIiDsupi9K5lJWVsX79egYNGlT9mNlsZtCgQaxatQqAvn37smXLFjIzMykqKuKbb75hyJAhNZ7zhRdeICgoqPojOjq6wb+P+iivtPDXzzbz8rdVG9uOu7AtU2917I1tG8vVPVrx50GxADz12WbdCSciImflskXpyJEjVFZWEh5++hYc4eHhZGdnA+Du7s6rr77KpZdeSkJCAo888sg573h78sknyc/Pr/7IyMho0O+hPgpLyrlzzjrmrs3AZILJQ+N4Zqjjb2zbmB78UywDO7WgtMLCvR+sJ/9kudGRRETEwbgbHcBow4YNY9iwYbU61svLCy8vx79kdSj/JGNnpbAjuxAfDzf+e3NPBjvhnm0NzWw28doNCVzz+s8cOHqCv8xP4+3bEjGZVCZFRKSKy44ohYaG4ubmRk5OzmmP5+Tk0LJlS4NSNbytWflcl1y1sW2ovxcfTeinknQOzfw8mXZrIp5uVXvC/boAp4iICLhwUfL09CQxMZHFixdXP2axWFi8eDFJSUkGJms4S3b+trFth7CqjW17RAUbHcvhdY8K4u/XdgXg5W93sGrPUYMTiYiIo3DqolRUVERqaiqpqakA7Nu3j9TUVNLT0wGYNGkS06dPZ86cOWzfvp17772X4uJixo4da2DqhvHhmnTumrOO4rJKktqF8Om9/Ylu7vwb2zaWm/pEV29z8sDcDWTnlxgdSUREHIBTz1Fat24dl156afXnkyZNAmDMmDHMnj2bG2+8kcOHD/PMM8+QnZ1NQkIC33777RkTvJ2ZxWLl5e92Mm3ZHgBG9IrkxRGuuWdbQzKZTPzftd3YmlXA9kMFTPxwA/Pu7oeHm/4eRUSasia911t9Gb3XW0l5JY/MT+OrTYcAeHhQRx68rIMmI9fD/iPFDH3jZwpLKhh3YVueGRpndCQREbEz7fXWBOQVlzF6xhq+2nQIDzcTr46K56FBsSpJ9RQT6sero+IBmLliH19uyjI4kYiIGElFyQntO1LMiDdXsP7AMQK83Zkzti/XJ7r+xraN5fKuLbl3YHsAHv9kE7/kFhqcSEREjKKiZIPk5GTi4uLo06dPo7/2uv15jHhzBfuPniAy2IcF9/anf4fQRs/h6h4Z3JGkdiEUl1Vyz/sbKC6tMDqSiIgYQHOU6qGx5yh9uSmLSR+nUVZhoUdUEDPG9CYswLvBX7epOlxYyjWv/0ROQSlD4yP4700JurQpIuICNEfJxVitVqYu3cP9H26krMLC4Lhw5t3dTyWpgbUI8OLN0b1wN5v4Ii2LOSv3Gx1JREQamYqSg6uotPDUwi289O0OAO7oH8O0WxPx9XTqlR2cRmKb5jx1dRcA/vHVdtYfOGZwIhERaUwqSg6sqLSCO+es48M16ZhM8Mw1cTw7rKs2tm1kd/SP4ZoeraiwWJn4wQaOFJUaHUlERBqJipKDys4vYdS0VSzbdRhvDzNv3ZrIuAFtjY7VJJlMJl66vgcdwvzJLijhwbkbqbRoap+ISFOgouSAtmUVcF3yCrYfKiDU35OP7k7i8q6uu5GvM/Dzcmfarb3w9XRj5Z6jvPr9TqMjiYhII1BRckDz12eQXVByamPbC4mPDjY6kgAdwgJ46foeALy5dA8/bMsxOJGIiDQ0FSUH9NRVXbj/0g58eo82tnU0Q+MjGHthDACTPk7lwNFiYwOJiEiDUlFyQO5uZv4ypBNBvh5GR5GzePLKLiS2aUZhSQX3vL+BkvJKoyOJiEgDUVGygZErc4vxPN3NJN/Si1B/T7YfKqheukFERFyPVuauh8ZemVscy7Jdhxkzcy0ebiYWTxpI6xBdJhURcQZamVukEVzSsQUXxYZSXmnltUW7jI4jIiINQEVJpB4eG9IZgIWpmWw/VGBwGhERsTcVJZF66B4VxNU9WmG1wivfaW0lERFXo6IkUk9/ubwTbmYTP+7IZe2+PKPjiIiIHakoidRT21A/buwTDcDL3+5A90eIiLgOFSURO3josli8PcysO3CMxdtzjY4jIiJ2oqIkYgfhgd6MvbBq0+JXvtupTXNFRFyEipKIndxzcXsCvd3ZmVPI/1IzjY4jIiJ2oKIkYidBvh7cd2kHAF79fhelFdraRETE2akoidjRmKQYwgO9yDx+kg/XpBsdR0RE6klFyQba601q4uPpxkOXdQTgjR9/oai0wuBEIiJSHypKNpg4cSLbtm0jJSXF6CjigG7oHUW7UD+OFpcx46e9RscREZF6UFESsTN3NzOPXN4JgOnL93K0qNTgRCIiYisVJZEGcGW3lnSPDKK4rJLkJXuMjiMiIjZSURJpAGazicevqNow9/3VBzh47ITBiURExBYqSiINZEBsKBd2CKGs0sJrP+w2Oo6IiNhARUmkAT02pGpUacHGg+zMLjQ4jYiI1JWKkkgDio8O5qruLbFaq7Y2ERER56KiJNLAHrm8E25mE4u257D+QJ7RcUREpA5UlEQaWPsW/tzQOwqAl77ZidWqDXNFRJyFipJII3jwsli83M2s3Z/H0p2HjY4jIiK1pKIk0ghaBflwR/8YAF76dgcWi0aVREScgYqSSCO5d2B7Arzd2ZFdyLyUDKPjiIhILago2UCb4ootgn09eeBPHQB49vOtrNuvid0iIo7OZNXMUpsVFBQQFBREfn4+gYGBRscRJ2CxWLnvgw18uzWbED9PFk68kOjmvkbHEhFpUury81sjSiKNyGw2MeXGeLpGBHK0uIy75qyjqLTC6FgiIlIDFSWRRubr6c6MMb1pEeDFzpxCHpq7kUpN7hYRcUgqSiIGaBXkw/Tbe+Plbmbxjlxe/naH0ZFEROQsVJREDJIQHcwro+IBeGv5Xj5epzvhREQcjYqSiIGGxUfw4GWxADz12WbW7tOdcCIijkRFScRgf74slqu7t6K80sqE99aRfvSE0ZFEROQUFSURg5nNJv41Kp7ukUEcO1HOnXNSKCwpNzqWiIigoiTiEHw83Zh+e2/CA73YnVvEA7oTTkTEIagoiTiIlkHeTL+9N94eZpbuPMw/v95udCQRkSZPRUnEgfSICuZfp+6Ee+fnfcxbm25wIhGRpk1FScTBXNMjgocHdQTgbwu3sGrPUYMTiYg0XSpKIg7owcs6MDQ+ggqLlbvfXceWzHyjI4mINEkqSiIOyGQy8crIHvSJaUZhaQW3vbOGndmFRscSEWlyVJREHJS3hxsz7+hDfFTVsgGjZ6xh7+Eio2OJiDQpKkoiDizA24M54/rSuWUAR4pKGT1jDRl5WpBSRKSxqCjZIDk5mbi4OPr06WN0FGkCgn09ef+uC2jfwo9D+SXcMmM12fklRscSEWkSTFarVava2aigoICgoCDy8/MJDAw0Oo64uOz8Em54axXpeSdo38KPjyYkEervZXQsERGnU5ef3xpREnESLYO8+eCuC4gI8mbP4WJunbGG4yfKjI4lIuLSVJREnEh0c18+GN+PFgFe7MguZMzMtdoXTkSkAakoiTiZtqF+fHDXBTTz9SDtYD7jZqdwoqzC6FgiIi5JRUnECXUMD+C9Oy8gwNudlP3HuPvd9ZSUVxodS0TE5agoiTipbpFBzB7bF19PN37+5QiPfrLJ6EgiIi5HRUnEiSW2acY7Y/rgbjbxRVoWX28+ZHQkERGXoqIk4uSS2odw78D2ADzzvy3kFetOOBERe1FREnEB9/+pA7Fh/hwpKuO5L7YaHUdExGWoKIm4AC93N14ZFY/ZBAtTs1i0LcfoSCIiLkFFScRFJEQHM/6idgA8tXAz+Se1vpKISH2pKIm4kIcHd6RdqB85BaU8/9U2o+OIiDg9FSURF+Lt4cbLI3tgMsHH6w6ybNdhoyOJiDg1FSURF9M7pjljkmIAePLTTdriRESkHlSURFzQY1d0Irq5D1n5Jbz07Q6j44iIOC0VJREX5OvpzksjegDw/up0Vu45YnAiERHnpKIk4qL6dwjllgtaA/DEp5u1ca6IiA1UlERc2JNXdiYiyJv0vBP867tdRscREXE6Kko2SE5OJi4ujj59+hgdReScArw9+OeI7gDMWrmPdfvzDE4kIuJcTFar1Wp0CGdVUFBAUFAQ+fn5BAYGGh1HpEZ/mZ/GJ+sP0i7Uj68fughvDzejI4mIGKYuP781oiTSBDx9dRxhAV7sPVLMa4t0CU5EpLZUlESagCBfD54fXnUJbvryvaRlHDc2kIiIk1BREmkiBseFMyw+AosVHv0kjdKKSqMjiYg4PBUlkSbk2WFdCfHzZFdOEck//mJ0HBERh6eiJNKENPfz5LlruwHw5tI9bM3KNziRiIhjU1ESaWKu6t6SK7q2pMJi5dH5myivtBgdSUTEYbnb8kWff/55nb9m8ODB+Pj42PJyImJHJpOJ567ryup9R9l2qIC3lu3h/j/FGh1LRMQh2bSOktlct4Eok8nE7t27adeuXV1fyqFpHSVxZp9tPMjDH6Xh6WbmywcH0DE8wOhIIiKNolHWUcrOzsZisdTqw9fX19aXEZEGcl1CJJd1DqOs0sKjn2yiQpfgRETOYFNRGjNmTJ0uo916660acRFxMCaTieeHdyfAy520jOPMXLHP6EgiIg5HW5jUgy69iSv4KCWdxz/djJe7ma8evIgOYf5GRxIRaVANeunt2LFj5OVVbax5+PBhFixYwNatW21LKiKGu6F3NBfFhlJaYWH0jNXsyC4wOpKIiMOoU1GaMWMGiYmJ9O7dm6lTpzJ8+HAWL17MTTfdxIwZMxoqo4g0IJPJxL9GxRMb5k9OQSmjpq5i1Z6jRscSEXEIdbr01qNHD9asWcPJkydp3bo1+/bto0WLFuTn53PJJZeQmpragFEdjy69iSs5fqKM8e+uI2X/MTzdzLx2YwJX92hldCwREbtrsEtv7u7u+Pj40Lx5czp06ECLFi0ACAoKwmQy2Z5YRAwX7OvJe3dewJCu4ZRVWrh/7gZmaYK3iDRxdSpKbm5ulJSUALBs2bLqx4uKiuybSkQM4e3hxpujE7mtXxusVvj7F9t44ZvtWCy650NEmqY6FaVFixbh5eUFVI0i/erEiRO8/fbb9k0mIoZwM5t47tquPDqkEwBvLdvLI/PTKKvQOksi0vTUaQuT35ej3wsLCyMsLMwugUTEeCaTiYmXdiAswIsnFmzms42ZBPl48OywrkZHExFpVHbZFLe8vJyMjAx27txZvXSAiDi/Ub2jef3mngB8lJJB/olygxOJiDQum4tSYWEhU6dO5ZJLLiEwMJCYmBi6dOlCixYtaNOmDePHjyclJcWeWUXEAFd2a0nnlgGcLK/k43UZRscREWlUNhWlKVOmEBMTw6xZsxg0aBALFy4kNTWVXbt2sWrVKiZPnkxFRQWXX345V1xxBbt377Z3bhFpJCaTiTv6xwDw7ur9VGpit4g0ITZtYXLzzTfzt7/9ja5dzz1fobS0lFmzZuHp6cm4ceNsDumotI6SNBUnyypJenExx0+UM/323gyOCzc6koiIzery81t7vdWDipI0JS98s523lu1lQIdQ3r/rAqPjiIjYrEH3ehORpum2fm0wm+DnX46wO6fQ6DgiIo2i3kXphRdeYObMmWc8PnPmTF566aX6nt4hJScnExcXR58+fYyOItJoopr5Vl9ym7Nqv7FhREQaSb2L0ltvvUXnzp3PeLxr165Mmzatvqd3SBMnTmTbtm26q0+anDGnJnV/uj6T/JNaKkBEGs53W7PZkV1gdIz6F6Xs7GxatTpz48wWLVpw6NCh+p5eRBxIUrsQOoVXLRUwX0sFiEgDWX/gGBPeW88V//7J6Cj1L0rR0dGsWLHijMdXrFhBREREfU8vIg7EZDJVjyq9u+qAlgoQkQaxNSvf6AjV6l2Uxo8fz5///GdmzZrFgQMHOHDgADNnzuThhx9m/Pjx9sgoIg7kup4RBPl4kJ53gqU7c42OIyIuqKLScX4Jq9Neb2fz6KOPcvToUe677z7KysqwWq34+Pjw+OOP8+STT9ojo4g4EF9Pd27sE83by/cye+V+LuuiNZVExL4cabS63iNKJpOJl156icOHD7Nq1SrS0tLIy8vjmWeesUc+EXFAvy4V8NPuI/ySq6UCRMS+KlypKAG888479OvXj4suuojevXuTmJjIjBkz7HFqEXFA0c19q0eS5qw8YHAaEXE1lRaL0RGq1bsoPfPMMzz00EMMHTqU+fPnM3/+fIYOHcrDDz+sUSURFzb216UCNhzkUP5JY8OIiEtxpBGles9Rmjp1KtOnT+fmm2+ufmzYsGH06NGDBx54gOeee66+LyEiDiipfQhdWgWy/VABw95YwdTRvegd09zoWCLiAlxqjlJ5eTm9e/c+4/HExEQqKirqe3oRcVAmk4m3bk2kc8sADheWcvP01XywRpfhRKT+TpZVGh2hWr2L0m233cbUqVPPePztt99m9OjR9T29iDiw1iG+fHpvf67q3pLySitPfbaFJxdsprTCcf6TExHnc7L8t/9DLAaPLtX70htUTeb+/vvv6devHwBr1qwhPT2d22+/nUmTJlUfN2XKFHu8nIg4ED8vd5Jv6cWbS/fwr+93MndtOjuzC5h2ayJhgd5GxxMRJ/T7olRptWLGZFiWehelLVu20KtXLwD27NkDQGhoKKGhoWzZsqX6OJPJuG9SRBqWyWRi4qUdiIsI5MG5G9mQfpxrXv+ZD8dfQIewAKPjiYiTKa347a63SosVDzfjstS7KC1ZssQeOUTEBVzaKYzP7x/A+HfX8UtuEQ/MTWXhxP54uRv4v5yIOB2r9bfLbUbfAWeXdZRERH7VNtSPD8dfQDNfD7YfKuDfi3YbHUlEnMzvepLhd8DZPKI0bty4Wh03c+ZMW19CRJxUWIA3L4zozj3vb+CtZXu4rHOYlg4QkVqz/K4pGV2UbB5Rmj17NkuWLOH48eMcO3asxg8RaZqu6NaKEb0isVhh0sdpFJVquRARqZ3flyOji5LNI0r33nsvc+fOZd++fYwdO5Zbb72V5s31G6OI/ObZYV1ZszeP9LwT/OPLbbx4fQ+jI4mIEyivdJyiZPOIUnJyMocOHeKxxx7jiy++IDo6mhtuuIHvvvvutElYItJ0BXp78K9R8ZhMMC8lg0XbcoyOJCJOoOJ3e71VGLzvW70mc3t5eXHzzTfzww8/sG3bNrp27cp9991HTEwMRUVF9sooIk4sqX0Id17YFoAnFmziaFGpwYlExNH9fkTJ6P1x7XbXm9lsxmQyYbVaqazUqrwi8pu/DOlEx3B/jhSV8dfPNhsdR0QcWEWlhbX78n773JlHlEpLS5k7dy6DBw+mY8eObN68mTfeeIP09HT8/f3tlVFEnJy3hxuv3ZiAu9nEd1tz2H6owOhIIuKg3l99+p6RFoOn89hclO677z5atWrFiy++yDXXXENGRgbz58/nqquuwmzW8kwicrquEUEMjgsH4ON1GQanERFHtTHj+Gmf//4ynBFsvutt2rRptG7dmnbt2rFs2TKWLVt21uMWLFhgczgRcS039I7mmy3ZLNyYyRNXdtaK3SJyBjfz6VueVThrUbr99tu1f5uI1MlFsaGEB3qRU1DK4u25XNW9ldGRRMTB/HEEqbTC2HnPNhel2bNn2zGGiDQF7m5mru8VxZtL9/DxugwVJRE5Q3nF6ZO3S8qdeDK3iEhd3dA7GoDluw5zKP+kwWlExNF0bhVw2udGjyjZrSitWbPGXqcSERcWE+pH37bNsVhhwYZMo+OIiIP545wklxlRGjVqlL1OJSIu7tdRpY/XZWAxeHsCEXEss1bsO+3zknInmqN0ww03nPVxq9VKXl7eWZ8TEfmjq7q35NnPt3Lg6AnW7s+jX7sQoyOJiIMoLju9GFUavI5SnYrSokWLeO+9985YTNJqtbJ8+XK7BhMR1+Xr6c7Q+FbMXZvBx+syVJRExGHVqSgNHDiQgIAALr744jOe69FDu4KLSO2N6h3N3LUZfL35EH8f1pUAbw+jI4mIIzL46nyd5igtWLDgrCUJ4IcffrBLIBFpGnpGB9MhzJ+ScgtfbjpkdBwRcTBe7o5xY369UmRnZ9srh4g0MSaTiRt6RwHa0kREqmQd/23JkM6tAgGwGjykVK+idPnll9srh4g0QcN7RuFuNrEx/Tg/7z5idBwRMdjXm38bXfZ2hRElq8Ez0Y2SnJxMXFwcffr0MTqKiFNrEeDFiF6RANz7/nq2ZRUYnEhEjPT7WhHg7X7GY0aoV1Fqqnu9TZw4kW3btpGSkmJ0FBGn99y13ejbtjmFpRXcMWstGXknjI4kIgb5dc2km/pEA47RMRxjXEtEmixvDzem396bTuEB5BaWMmbWWvKKy4yOJSIGOHGqKPl4ulU/ZvS1KxUlETFckI8Hc8b1JSLIm72Hixk3O4UTZRVGxxKRRnby1GKTvp5uOMpFq3oVJTc3t/MfJCJSCy2DvHn3zr4E+XiQmnGcBz7cqO1NRJqY2Sv3A1WL0jqKehWljRs32iuHiAgdwgKYeUdvvNzNLN6Ry9yUdKMjiUgj+f0NYrtyCn/3uBFpfqNLbyLiUBLbNOeJKzsD8PK3OzlSVGpwIhFpDAUlv11udzObHGQqt4qSiDig2/q1Ia5VIPkny3nxmx1GxxGRRnD0d78UeZh/qydOveAkQEpKCpdddhk9evRgxIgRPPfcc3z++eekp2vIXERs4+5m5h/Du2EywSfrD7J2X57RkUSkgR0p+u1u1yu6t3SNydwAt912G25ubtx99920bduWZcuWMXbsWGJiYggJ0Y7gImKbXq2bcVOf1gD8beFmyistBicSkYaUW1gCgLvZxMCOLaofN3qOUr2nlWdkZPDVV1/Rvn370x4/cOAAqamp9T29iDRhjw3pxHdbs9mVU8SsFfu4++L25/8iEXFKB45WLTY7LD4Ck8mEo8xSqveI0oUXXsjBgwfPeLxNmzZce+219T29iDRhzfw8qyd2/3vR7tM2zBQR15J/shyA0ACv0x43epEQm4rSiBEjePbZZ/nss8+45557+L//+z+OHTtm72wiIozsFUXvNs04UVbJwx+lUlyqhShFXNGvi016e1St0ejUc5Tat2/PihUrmDBhAiNHjuTHH3+kY8eO3HXXXcyYMYP169dTVqYtCESk/sxmE88P746vpxtr9uVxy/TV2uJExAWd/HX7Eg/HWszapqL0yiuv8MMPP5Cbm0tGRgZffPEFf/7zn8nPz+ell16ib9++BAQE0KNHD3vnFZEmqFPLAD4c349gXw/SDuYzatpKMnUZTsSl/FaU/lBNDJ7NXe/J3JGRkURGRnL11VdXP1ZUVERqaippaWn1Pb2ICAAJ0cF8ck8St72zlj2Hixk5dSXv3dmXDmEBRkcTETs4mFc1mbuZnyfg5JfezrdGkr+/PwMGDGDixIkAZGZm2vIyIiKn6RAWwCf39qddCz8O5Zdww1urq28pFhHnlpVf9W+5Q5j/aY875WTuPn36MGHCBFJSUmo8Jj8/n+nTp9OtWzc+/fRTmwOKiPxeZLAPn9zTn07hAeQVl/Hqd7uMjiQidnDi1I0a/l5VF7scZXkAmy69bdu2jeeff57Bgwfj7e1NYmIiEREReHt7c+zYMbZt28bWrVvp1asXL7/8MldddZW9c4tIE9bcz5N/jujG9VNX8fH6DG5LakO3yCCjY4mIjY4WlVJ86q43X8/Tq4nRC07aNKIUEhLClClTOHToEG+88QaxsbEcOXKE3bt3AzB69GjWr1/PqlWrVJJEpEEktmnO0PgIrFb4vy+3nbbzuIg4l3kpGdV/9vM6ddebYwwo1W8yt4+PDyNHjmTkyJH2yiMiUmuPX9GJ77dms2ZfHt9tzeGKbi2NjiQiNvDz/G1JgD8uD2D0L0H1XplbRMQoUc18GX9ROwD++fV2SisqDU4kIrY4empttJv7tsZ06nY3BxlQsq0obdq0CYul9htUbt26lYoKraYrIvZ378D2tAjwIj3vBLNX7Dc6jojYIDXjOACtgrzPeM7oi+o2FaWePXty9OjRWh+flJR03iUFRERs4eflzqNDOgHw+o+/sPdwkcGJRKSu0k+toRTXKrD6MZODLKRk0xwlq9XK008/ja+vb62O13YmItKQRvaK4r1VB9icmc+1b6zgXzfEM6Sr5iuJOItfN8RtE1K7XtGYbCpKF198MTt37qz18UlJSfj4+NjyUiIi52U2m5gxpjcTP9jAugPHmPDeeiZc0o5HL++Eu5umYoo4MqvVSsGpohTk43GW5xs70elsKkpLly61cwwRkfoJD/Rm7t39ePGbHbzz8z7eWraX1PTjvH5LT8ICzpz3ICKOoai0AsupMhT4u6LkGBfedNebiLgQDzczT18Tx5uje+Hn6caafXlc89+fSdmfZ3Q0EalBQUnVzV6e7ma8/7A0ADjpZG4REUd2VfdWfP7AAGLD/MktLOWmt1cz46e9hq/HIiJnOnZqaYBA79MvuznIXG4VJRFxTe1b+LNw4oUMi4+g0mLlH19t55H5aSpLIg5m/9FiAKKbn30us9H/ZlWURMRl+Xm585+bEnju2q64m00s2JDJyj21X9pERBre7pyqJT1iw/xPe9xBBpRUlETEtZlMJm5PiuHWfm0AePX7nYb/hioivzl47CQAbUL8DE5ydnYpSuXl5WRkZLBz507y8jRpUkQcz30D2+PlbmZD+nGW7TpsdBwROeVocSkALfy9TnvcURactLkoFRYWMnXqVC655BICAwOJiYmhS5cutGjRgjZt2jB+/HhSUlLsmVVExGZhgd7cnlQ1qvTaD7s0qiTiII4WVU3mDvH3NDjJ2dlUlKZMmUJMTAyzZs1i0KBBLFy4kNTUVHbt2sWqVauYPHkyFRUVXH755VxxxRXs3r3b3rlFROpswiXt8fFwI+1gPou35xodR6TJs1qtbM7MByD0DyNKvx3TmInOZNOCkykpKSxfvpyuXbue9fm+ffsybtw4pk2bxqxZs/jpp5+IjY2tV1ARkfoK9fdiTP8Ypi3bw5QfdnFZlzCHGd4XaYp25hRW/zmq2el3vTnKv0ybitLcuXOr/5ybm0tYWNhZj/Py8uKee+6xLZmISAOYcHE73lu1n22HCvh6czZX92hldCSRJmv/kRPVfw6paUTJ4CUn6z2Ze+TIkVRWVp71uYqKivqeXkTErpr5eXLnRe0AmPz5Vo4WlRqcSKTpOnisqiid9RcWBxlSqndRCg4O5sEHHzzj8aNHjzJo0KD6nl5ExO7uG9ie2DB/jhSV8tfPNmtit4hBNmYcB6B9aM1LAxj9z7PeRendd9/lhx9+YObMmdWPbd++nb59++Ln55hrIohI0+bt4cZrNybg4Wbiu605fLoh0+hIIk1SRl7ViFL3qOAznjM5yJCSXUaUPv30Ux599FHWrl3Ld999R1JSEtdddx1ffPGFPTKKiNhdt8gg/jyoIwDPfr61+j9sEWk8mw5W3fHW3M+jxmOMHu+1aTL3iBEjSEhIqP7o3r07b7zxBldddRUlJSW8/vrrjB071t5ZRUTsasLF7fhxRy7rDxzjL/PTmDu+H2azY/wWK+LqDuWfrP5zq6Az93lzlBtSbRpRat++PT/99BN33XUXMTExhISEMH36dKxWK7fccgu9evWivLzc3llFROzK3c3MlBvi8fFwY82+PD7ZcNDoSCJNRnZ+SfWfI4LPviGuI7BpROmVV16p/nNmZiapqamkpqYSEhLCkiVLeOedd3B3d6dz586kpaXZLayIiL21CfHj4cGx/PPrHbzw9XYGdQmnuZ9jrhAs4kqOnahakbtbZOA5jzN6MrdNRen3IiMjiYyM5Oqrr65+rKioiNTUVJUkEXEKYy9sy4INmezILuTFb7bz8sh4oyOJuLxfty5p7nf29ZMc5MqbbZfe0tPTz/m8v78/AwYMYOLEiUDVqJOIiKPycDPz/PBuAHy87iBr92lzb5GG9vinmwDYf6T4nMc55YKTffr0YcKECefc9DY/P5/p06fTrVs3Pv30U5sDiog0hsQ2zbm5bzQAT322mbIKi8GJRFyb5VT/aRFQw4iSgwwp2XTpbdu2bTz//PMMHjwYb29vEhMTiYiIwNvbm2PHjrFt2za2bt1Kr169ePnll7nqqqvsnVtExO4ev6Iz323NYXduEU8s2MSro+K1F5xIA2nm68GxE+U8dXWXcx5n9Bwlm0aUQkJCmDJlCocOHSI5OZnY2FiOHDnC7t27ARg9ejTr169n1apVKkki4jSCfT157cYE3MwmFmzI5PUffzE6kohLKimv5NiJqrvj24f6n/UYR1lwsl6TuXNzc/H09OSWW26hb9++9sokImKYSzq24Llru/LUZ1uY8sMuWjf35bqekUbHEnEpaae2Lgn19yTQp973lTUom1fmnjt3Lh07duTaa68lKSmJ3r17c/jwYXtmExExxOgL2jDh4qqNcx/7ZJMmd4vY2Y7sQgC6tAqs8fK2o1z1trko/f3vf+eWW25hx44dfP/99wA88cQTdgsmImKkx6/ozJXdWlJWaeHu99ax93CR0ZFEXMbzX28HoKCk4rzHGr1ptc1Fae/evUyePJmOHTty2WWX8f777zNv3jx7ZhMRMYzZbGLKDQnERwdz/EQ542ankFdcZnQsEZfg6+kGQJvmvjUe4/QjShUVFfj6/vYNdu7cGYvFQnZ2tl2CiYgYzcfTjRm39yaqmQ/7j57g7nfXUVJeaXQsEacX6F21Ce4tF7Q2OMn52VyUAObMmcPKlSspKqoaknZ3d+fECe3ALSKuo0WAF7Pu6EOAtzvrDhzjtUW7jI4k4tRyC0pIz6vqCq3PMaL0K6dcHgDgoosu4h//+AcDBgwgODiY2NhYSkpKeOedd1iyZAmFhYX2zNlghg8fTrNmzRg5cqTRUUTEQcWGBzDlhgQA3vlpHzuzneP/NxFHtP3Uv592LfzOsxmuY1x7s7koLVu2jPz8fHbu3Mn777/P8OHDueSSS5g6dSqXXXYZzZo1o0uXcy8i5Qgeeugh3n33XaNjiIiDGxwXzuVx4VRYrPxt4WYsFoN/zRVxUntyq65CxYadff2kPzL6X1q9Fy+IjY0lNjaWm266qfqxffv2sW7dOjZu3Fjf0ze4gQMHsnTpUqNjiIgTmDysKz//coSU/cf4ZMNBbugdbXQkEaez90hVUWrX4txFyeknc59L27ZtGTVqFP/85z/rdZ7ly5czdOhQIiIiMJlMLFy48IxjkpOTiYmJwdvbmwsuuIC1a9fW6zVFRGoSGezDnwfFAvDC19s5WlRqcCIR55N6arHJjuG1HFFy1jlKjaG4uJj4+HiSk5PP+vxHH33EpEmTmDx5Mhs2bCA+Pp4hQ4aQm5tbfUxCQgLdunU74yMrK6uxvg0RcSFjL2xL55YBHDtRzoT31usuOJE6KK+0sONQ1RylC9qGnPNYBxlQqv+lt4Z05ZVXcuWVV9b4/JQpUxg/fjxjx44FYNq0aXz11VfMnDmzevHL1NRUu+UpLS2ltPS33yALCgrsdm4RcQ4ebmbeuKUnI95cyboDx/jL/DT+e1NPzGZH+W9dxHGtP3CMCouVAG93WgZ61+prrAbPUnLoEaVzKSsrY/369QwaNKj6MbPZzKBBg1i1alWDvOYLL7xAUFBQ9Ud0tOYniDRFHcICmHZbIh5uJr7cdIjJn29ld06hJniLnMf9H1bNXfZyN5/3lwuXnqPUGI4cOUJlZSXh4eGnPR4eHl6nRS8HDRrEqFGj+Prrr4mKijpnyXryySfJz8+v/sjIyLA5v4g4t/7tQ3lxRA8A3lt9gMGvLafH378neckvBicTcVxHTs3rO3ai3OAktefQl94aw6JFi2p9rJeXF15eXg2YRkScyfWJUQB8tC6DzQfzKSqt4L+LdzOmfwz+Xk3+v1eR0/x+z7b37uxbh69riDS157QjSqGhobi5uZGTk3Pa4zk5ObRs2dKgVCLS1FyfGMXHE5LY/OzltAv1o7TCwuLtOef/QpEm5ujv9kpMbNPsvMebHGQ6t9MWJU9PTxITE1m8eHH1YxaLhcWLF5OUlGRgMhFpitzdzFzdoxUAX6QdMjiNiOPJOn4SgLAAL7zc3Wr9dUbP/HPoolRUVERqamr1nWv79u0jNTWV9PR0ACZNmsT06dOZM2cO27dv595776W4uLj6LjgRkcZ0TY8IAJbvOkxBifPMwRBpDGkH8wGIanaubUt+4yiTuR36Ivq6deu49NJLqz+fNGkSAGPGjGH27NnceOONHD58mGeeeYbs7GwSEhL49ttvz5jgLSLSGDqG+9MhzJ9fcov4YWtO9RwmEYGZP+8Dqm6EqBODJyk5dFEaOHDgaZO/zub+++/n/vvvb6REIiI1M5lMXNOjFf9etJsvN2WpKImcUlFpYd+RYgBia7kit4MMKDn2pTcREWdzzal5Sj/tPkK+E90CLdKQduYUVv/5ym6t6vS1mqPkhJKTk4mLi6NPnz5GRxERB9MhLIDOLQOosFh58dvtVFRajI4kYrhf93cb0CEUT/faVQ+Tg0xSUlGywcSJE9m2bRspKSlGRxERB3TvwPYAzF2bwV3vrqNQE7uliUvZlwdAQnSwsUFsoKIkImJn1yZEMnV0L7w9zCzdeZib3l6tsiRNVlmFhe+3Va0tlhhz/vWT/kgLToqIuKAru7fio7uTCPX3ZGtWAfe+v4GyCl2Gk6ZnV04hJ8oq8fN045LYFkbHqTMVJRGRBhIfHczssX3x9XTj51+O8MSCTee9k1fE1bz2wy4AEloHn3cj3LOxGjydW0VJRKQBdYsMInl0L9zMJhZsyOS5L7epLEmTUVhSzuIduQBcUce73RxkLreKkohIQ7u0UxgvDO8OwKwV+/n7FypL0jSs23+s+s/D4iNsOofR/1RUlEREGsENfaJ5cUR3TCaYvXI/976/gZ3Zhef/QhEn9tK3OwC4tV9rgnw86vS12hRXRKSJualva166vgcmE3y7NZsh/17Og3M3Uq61lsQFfZySwY5TvwwM6drS5vMYPfaqoiQi0ohu6B3NF/cP4OrurTCZ4PO0LD7bkGl0LBG7qrRYeezTTdWf94lpXudzaI6SiEgT9esE7yeu6AzA1GV7qLQY/XuziP2s2Xu0+s+v3RiPt4ebzefSHCUnpC1MRMQebu3XhmBfD/YdKearzYeMjiNiN99syQYgMtiH4T1t2xzaQQaUVJRsoS1MRMQe/LzcGdu/LQBvLvkFi0aVxAVYLFY+21h1OflvV3cxOE39qSiJiBjojv4x+Hu5syO7kLkp6UbHEam3jRnHKCqtwN/LncFx4fU+nxacFBFpwoJ8Pao30X164Ra+25ptcCKR+lm8vWqByUs6tsDdzfaaocncIiICwH0D23ND7ygsVnhg7kZW/24irIgzsVqtvLl0DwCJbeq+Ae7ZT2qf09hKRUlExGAmk4l/Du/O5XHhlFVYuOf99ew7Umx0LJE6O36ivPrPg7rU77KbyUGGlFSUREQcgLubmf/e3JP46GCOnyjnztkpHD9RZnQskTr5X2rVJO7mfp60DvG1yzmNvsVBRUlExEF4e7gx/fZEIoK82XukmL/M36Q94cSp7D01Ehoe6F3vcznGeJKKkoiIQwkL8Gb6mN54uJlYtD2nej0aEWfw60TuB//UwW7nNPqXBRUlEREH0zUiiHsvqboTbvLnW8n/3bwPEUeVW1BC5vGTmE1wcccW9T+hgwwpqSiJiDig+y7tQLsWfhwuLGXoGz/z1Geb2Z1TaHQskRqtP3AMgNiwAPy83A1OYz8qSiIiDsjbw42Xr++Bn6cb6Xkn+GBNOje+vZqDx04YHU3krL48tQ3PxR1D7Xpeo6fpqSjZQHu9iUhj6B3TnBVP/Im3bkukS6tA8orLmPDeek6WVRodTeQ0JeWVLN6eA8Cw+Ei7nNPkINfeVJRsoL3eRKSxBPt6MqRrS2aM6U1zP0+2ZhUw4f31WjpAHMqG9GOUlFsIC/CiW2SgXc9t9H2fKkoiIk4gMtiHN0f3wsvdzPJdh7nm9Z81Z0kcxqo9VavJJ7UPsdtCkQ6y3qSKkoiIs+jXLoQF9/WndXNfDh47yVMLtxgdSQSAJTurlgXo3z7Ebufs1boZY5La0CfGTluh2EhFSUTEiXSNCOKjCf1wN5tYuy+PbVkFRkeSJm5XTiFbMqvehxd2sN9E7sFx4fz92m5c0a2V3c5pCxUlEREn0yrIhyHdWgIwZ+V+Y8NIk5eWcRyAnq2DiWpmn21LHImKkoiIE7qjfwwAC1Mz2XwwX5O7xTA//3IEgKR29rvs5khUlEREnFDvNs3oGhFIaYWFoW/8TJ/nFzFvbbrRsaSJqbRYWbbrMACXdg4zOE3DUFESEXFCJpOJx6/oTOvmvgR6u1NeaeWJBZt5d9V+o6NJE5J28DjHT5QT6O1Oz+hgo+M0CBUlEREndXHHFix/7FLSJl/OXQPaAvDM/7by8+4jBieTpmLpjqq73S7q2AJ3N9esFK75XYmINCEmk4mnru7CzX2jAXjskzQKSrSRrjS8pacuuw20xya4DkpFSUTEBZhMJv52dRxtQnzJyi/hyQWbqai0GB1LXNjhwlI2HcwH4JJOKkoiIuLg/Lzc+deoeMwm+GrTIe6cs46i0gqjY4mLWn5qNKlbZCBhAd4Gp2k4Kko20Ka4IuKo+sQ0563beuPtYWbZrsM8/FEqVqO3XxeX9HlaFgCXdnLNu91+paJkA22KKyKObHBcOB/cdQEebiZ+2JbDvJQMoyOJizmUf7J6WYCruhu7cnZDU1ESEXFBiW2a8+iQTgBM/nwrD83byNasfINTiatYtC0HqFqNu0urQIPTNCwVJRERF3XXgHYMjgunrMLC/1KzuOnt1eQWlhgdS1zAV5sPAXDlqa10XJmKkoiIizKbTbx1ayKf3tufLq0CKSyp4B9fbjc6lji5w4WlrN2XB8CVBm9Y2xhUlEREXJjZbCKxTTNevr4HZlPVBNzXftilu+HEZv9LzcRihfioIKKbu94muH+koiQi0gR0jwrirovaAfCfxbu57NWlrD9wzOBU4ox+vTlgZO9og5M0DhUlEZEm4okrOvOfmxJoE+JLTkEpN729ivdW7dfyAVJrO7ML+SW3CA83E8PiI4yO0yhUlEREmgiz2cS1CZF8/eBFXN29FeWVVp7+31Ye+2QTFovKkpzf52mZAFzSMYwgHw+D0zQOFSURkSbGz8udN27pyVNXdcHNbGL++oN8suGg0bHEwVksVhZurFpk8tqEpjGaBCpKIiJNkslkYvzF7Xjs1FpLL36zg+MnygxOJY5s9b6jZB4/SYCXO4Pjwo2O02hUlEREmrBxA9oSG+ZPXnEZd85ZR3a+1lmSs/tkXdWo4zXxrfD2cDM4TeNRURIRacI83My8PLIHAV7urD9wjGFv/Myew0VGxxIH80tuEV+eWmTyhiZyt9uvVJRERJq4nq2b8cUDA+gY7k9uYSk3v72avSpL8jtzVu6nrMLCRbGhJEQHGx2nUakoiYgIMaF+zLs7ic4tA6rK0vTV7DtSbHQscQB5xWV8emqy/x39YzCZTAYnalwqSiIiAkBzP08+uOsCOoUHkFNQyqhpK7UopfD28r2cKKukW2Qgl3YKMzpOo1NRskFycjJxcXH06dPH6CgiInYV4u/FB+MvIK5VIEeKyhg1bSW3vbOGXTmFRkcTAxSXVjAvJR2AB/8Ui9nctEaTAExWLclqs4KCAoKCgsjPzycwMNDoOCIidlNcWsFjn27iq01VE3jDA734/P4BhAd6G5xMGtPfv9jKrBX7iQz2YdmjA3F3c43xlbr8/HaN71hEROzKz8ud5Ft6sfQvA4kN8yenoJSr//szr36/k4pKi9HxpBHknyxn3tqqfd3+Mbyby5Skumqa37WIiNRKTKgf74zpQ6sgb44UlfL6j78wf71W8W4Kpi7dw8nySjqG+zOwYwuj4xhGRUlERM6pdYgvSx8dyAN/6gDAm0t/oVyjSi5t1Z6jTFu2B4BHh3Rucne6/Z6KkoiInJeXuxv3DexAiJ8nGXknufDFH5mzcr/RsaSB/HfxbgASooMZ1KXp3en2eypKIiJSKz6ebjw0KBaA3MJSnv1iK+sP5BmcSuxt6tI9rNp7FIBXb4hv0qNJoKIkIiJ1cHtSDGv/ehnDe0ZitcLEDzby3qr9FJdWGB1N7CC3sIR/L9oFwDU9WtG+hb/BiYynoiQiInUSFujNs0O70rq5L9kFJTz9v60kvbCYf369nbziMqPjiY2sVisT3ltPaYUFTzczzwyNMzqSQ1BREhGROgvy9eDrhy7i78O6EhPiS0FJBW8v38utM9ZwsqzS6Hhig5kr9rMx/TgAM8b0JixAa2aBipKIiNjI38udMf1j+PGRgbwzpjchfp5sO1TAUws3Gx1N6mjP4SKmfL8TgAf+1IGLm/ByAH+koiQiIvViNpu4rEs4yaN7YTbBgg2Z/Lgjx+hYUks7swu58a3VFJdVkhAdzMODOhodyaGoKImIiF30axfCuAvbAjBu9jpueGsVRZrk7dCsVisj3lzBkaJSIoN9mH577ya5n9u5qCiJiIjdTLq8I/3bhwCwdl8eM37aS6VFW4o6IovFyp1z1lF8ak7ZjDG9aRHgZXAqx6OiJCIiduPr6c6H4/vx6qh4AP69aDddnvmWjenHDE4mfzTlh138uCMXgImXtqdLK23ufjYqSiIiYnfDe0aS2KYZAGUVFiZ/vlV3wzmIkvJK/jxvI28s+QWAV0fF8+iQzganclwqSiIiYndms4n377yAD++6ADeziU0H8+nyzLfMW5tudLQmb8oPu1iYmgXAiJ6RjOgVaXAix6aiJCIiDcLH043+HUL5y+Wdqh+b/PlWvt+arXlLBvlq0yHeXr4XgDv6x2iLklpQURIRkQZ178D2bH/uCvq3D6G0wsLd761n9IzVpGUcx2pVYWosFZUWXv5uBwATLmnHs8O6qiTVgoqSiIg0OB9PN94c3Ys7+sfg6+nG6r15XJu8ggfnpWLR6FKj+GFbDgeOnqCZrwcP/inW6DhOQ0XJBsnJycTFxdGnTx+jo4iIOI1gX0+eHdaVT+/tz0WxoQB8kZbF2NkpTFu2h21ZBQYndG0zV+wDYPQFbfDzcjc4jfMwWTXuabOCggKCgoLIz88nMFC3VYqI1MXCjZk8/HEqv/4UMpng78O6cntSjKG5XNH6A3lcP3UVnm5mfnr8UsIDm/Y+bnX5+a0RJRERMcR1PSOZO74ff+ocRrfIQKxW+MdX29lzuMjoaC5n6tKqCdzDe0Y2+ZJUVypKIiJimH7tQph5Rx++uH8AF3dsQVmFhbGzUjh47ITR0VxGyv48Fm3PwWSC8Re3MzqO01FREhERw5lMJl4Z2YPWzX1JzzvBNa//zNy16VpGoJ4qLVb+/sVWAG7qE02HMH+DEzkfFSUREXEI4YHezLu7Hz2igjh+opwnF2zmhrdW6VJcPXy2MZMtmQUEeLnzyO/Ws5LaU1ESERGHERHsw/x7knj6mjj8vdxZf+AYV//3J1b+csToaE7HarXy3qr9ANwzsD2h/trw1hYqSiIi4lC83N24c0Bbvnv4Yvq1a05JuYXbZq7lL/PTyMjT3KXa+jwti7SD+Xh7mBnVO8roOE5LRUlERBxSZLAPs8f25erurai0WPlk/UEGTVnG28v3sGbvUa3qfQ75J8v5vy+3AzBxYAfCAnSnm61UlERExGF5e7iRPLoXC+7rT1K7qi1Q/vn1Dm58ezVP/28LpRWVWtn7LP762WaOFJXSvoUfd1+iO93qQ0VJREQcXq/Wzfhw/AWMu7Bt9WPvr06n09++pd1fv+Y/i3YbmM6xbM3K56tNhzCZ4OWR8Xi5uxkdyampKImIiFMwmUw8MzSO1U9exuNXdMbd/NuGrq8t2sWunEID0zmGo0WljJ2VAsDgLuEktmlmcCLnp6IkIiJOpWWQN/cObM/SRwfSLfK37Scuf205q/YcNTCZsaxWK//8ege5haVENfPhmaFxRkdyCdrrrR6015uIiPG+35rN3e+tr/78qu4tubRTGAHeHlzRraWByRrXlB928d/FVZcgp92a2KS+97qqy89vFaV6UFESEXEM+SfKGfzaMnILS097fOroXlzZvZVBqRrP0p25jJ2dgtUK9w5sz2NDOmEymc7/hU2UilIjUVESEXEcx4rLSM04zsMfp3L8RDkAnu5merUOJibEj39c1w13N9ebcZJ5/CRDXltOUWkFw3tG8tqNCUZHcnh1+fnt3kiZREREGlQzP08u7RzGrDv6sGBDJluy8tmYfpzVe/NYvTcPbw83/jwolmBfT6Oj2s2O7ALGzkqhqLSC2DB/XhjR3ehILkcjSvWgESUREcdlsVj5+ZcjvPHjL6zdnweAt4eZyzqHMzgunM6tAujc0jn/7660WHnn5738d/EvFJVW0L6FH7PH9iW6ua/R0ZyCLr01EhUlERHHZ7Va+XRDJjN+2suO7NOXEBiT1IYHLot1un3Qnlywmblr0wHo27Y5b9+W6FIjZQ1NRamRqCiJiDgPq9XKqr1Hue+DDdVzmACCfT3o3aYZ/dqFcOeAtg49CbqkvJJ/fr2dd1cdAGBkYhTPD++mRSXrSEWpkagoiYg4n905haTsP0aIvycvf7uDPYeLT3v+xt7RvDSyh0Hpzs5qtfL8V9v5YE06J8srAXjyys5MuKS9wcmck4pSI1FREhFxbrmFJTzx6WZ+3JF7xnMvXd+dQG8PktqHGHpZa+nOXO44tdo2QJCPB/93XTeGxUcYlsnZqSg1EhUlERHXseKXI4yeseaMxwO83BkQG0qfmObc0T8Gs7nhL82dKKtg5s/7WLAxk72/G/Eak9SGp66Ow9Pd9ZY5aEwqSo1ERUlExLWcKKtg4cYsth3KZ/PBfI4Wl3Hw2Mnq50P9Pbk9KYb2Lfw5WlzKviPFjLuwrd3uNjtwtJh5KRm8v+oAhaUVpz338YQk+rZtbpfXaepUlBqJipKIiGurqLQwLyWDKT/sIq+4rMbjhnQNx91sZlTvKBLbNCPA2+O857ZarVRYrMxbm87y3UfYdPA4OQWnryx+XUIEj17Rmchgn3p/L/IbFaVGoqIkItJ0HC0q5ZP1B9mSVcCavUfP2C7lV+5mE51aBrA1q4DWzX0J9vVgS2Y+lt/9tO0RFcT2QwWUV57+I9hkgqhmPvRp05yHB3fUukgNREWpgSUnJ5OcnExlZSW7du1SURIRaaK+3ZJNyv48luzMPW0uUV1dHhdOfHQwIxOjCA/0tmNCORsVpUaiESUREYGqy2h5xWUUl1ayO7eQLZkFrN1/lJ7RzWjm58lHKem4m80UlVaQnneCG3tHkxjTjIToYIJ9PQgLUDlqTCpKjURFSURExPnU5ee37i8UERERqYGKkoiIiEgNVJREREREaqCiJCIiIlIDFSURERGRGqgoiYiIiNRARUlERESkBipKIiIiIjVQURIRERGpgYqSiIiISA1UlERERERqoKIkIiIiUgMVJREREZEaqCiJiIiI1MDd6ADOzGq1AlBQUGBwEhEREamtX39u//pz/FxUlOqhsLAQgOjoaIOTiIiISF0VFhYSFBR0zmNM1trUKTkri8VCVlYWAQEBmEwmAPr06UNKSsp5v7Y2x53vmIKCAqKjo8nIyCAwMLBu4R1cbf8enfH17XFuW89R16+ry/H1fU+78vsZXPc9ba/z6j3tfIx8T9f3ta1WK4WFhURERGA2n3sWkkaU6sFsNhMVFXXaY25ubrX6B1Gb42p7rsDAQJf7R1jb790ZX98e57b1HHX9urocb6/3tCu+n8F139P2Oq/e087HyPe0PV77fCNJv9JkbjubOHGi3Y6r7blckdHfe0O+vj3Obes56vp1dTle7+lzM/p7b6jXt9d59Z52PkZ+74352rr05sQKCgoICgoiPz/fJX9bkaZF72dxNXpPuwaNKDkxLy8vJk+ejJeXl9FRROpN72dxNXpPuwaNKImIiIjUQCNKIiIiIjVQURIRERGpgYqSiIiISA1UlERERERqoKLUBGRkZDBw4EDi4uLo0aMH8+fPNzqSSL0NHz6cZs2aMXLkSKOjiNTZl19+SadOnYiNjWXGjBlGx5Fz0F1vTcChQ4fIyckhISGB7OxsEhMT2bVrF35+fkZHE7HZ0qVLKSwsZM6cOXzyySdGxxGptYqKCuLi4liyZAlBQUEkJiaycuVKQkJCjI4mZ6ERpSagVatWJCQkANCyZUtCQ0PJy8szNpRIPQ0cOJCAgACjY4jU2dq1a+natSuRkZH4+/tz5ZVX8v333xsdS2qgouQAli9fztChQ4mIiMBkMrFw4cIzjklOTiYmJgZvb28uuOAC1q5da9NrrV+/nsrKSqKjo+uZWqRmjfmeFmls9X1/Z2VlERkZWf15ZGQkmZmZjRFdbKCi5ACKi4uJj48nOTn5rM9/9NFHTJo0icmTJ7Nhwwbi4+MZMmQIubm51cckJCTQrVu3Mz6ysrKqj8nLy+P222/n7bffbvDvSZq2xnpPixjBHu9vcSJWcSiA9bPPPjvtsb59+1onTpxY/XllZaU1IiLC+sILL9T6vCUlJdaLLrrI+u6779orqkitNNR72mq1WpcsWWK9/vrr7RFTxCa2vL9XrFhhve6666qff+ihh6wffPBBo+SVutOIkoMrKytj/fr1DBo0qPoxs9nMoEGDWLVqVa3OYbVaueOOO/jTn/7Ebbfd1lBRRWrFHu9pEUdVm/d337592bJlC5mZmRQVFfHNN98wZMgQoyLLeagoObgjR45QWVlJeHj4aY+Hh4eTnZ1dq3OsWLGCjz76iIULF5KQkEBCQgKbN29uiLgi52WP9zTAoEGDGDVqFF9//TVRUVEqWeIQavP+dnd359VXX+XSSy8lISGBRx55RHe8OTB3owNIwxswYAAWi8XoGCJ2tWjRIqMjiNhs2LBhDBs2zOgYUgsaUXJwoaGhuLm5kZOTc9rjOTk5tGzZ0qBUIrbTe1pcmd7frkdFycF5enqSmJjI4sWLqx+zWCwsXryYpKQkA5OJ2EbvaXFlen+7Hl16cwBFRUX88ssv1Z/v27eP1NRUmjdvTuvWrZk0aRJjxoyhd+/e9O3bl3//+98UFxczduxYA1OL1EzvaXFlen83MUbfdidVtzgDZ3yMGTOm+pjXX3/d2rp1a6unp6e1b9++1tWrVxsXWOQ89J4WV6b3d9Oivd5EREREaqA5SiIiIiI1UFESERERqYGKkoiIiEgNVJREREREaqCiJCIiIlIDFSURERGRGqgoiYiIiNRARUlERESkBipKIiIiIjVQURIRERGpgYqSiIgDGD58OM2aNWPkyJFGRxGR31FREhFxAA899BDvvvuu0TFE5A9UlETEofzlL3/huuuuq9WxAwcOxGQyYTKZSE1NtekcjmLgwIEEBASc9bk77rij+vtcuHBh4wYTaeJUlETEoaSmppKQkFDr48ePH8+hQ4fo1q3baefo0aNHjV/za/G45557znhu4sSJmEwm7rjjjrrEblD/+c9/OHTokNExRJokFSURcShpaWl1Kkq+vr60bNkSd3f3084RHx9/zq+Ljo5m3rx5nDx5svqxkpISPvzwQ1q3bl3n3OeTkJBAt27dzvjIyso679cGBQXRsmVLu2cSkfNTURIRh3Hw4EGOHDlSXZSOHz/O0KFDGTBgANnZ2XU6B8DgwYPx9fWlU6dOrFmz5rTjevXqRXR0NAsWLKh+bMGCBbRu3ZqePXueduzAgQO5//77uf/++wkKCiI0NJSnn34aq9VafYzFYuHll1+mQ4cOeHl50bp1a55//vnq51NTU9myZcsZHxEREXX6OxKRxqWiJCIOIzU1leDgYGJiYti8eTN9+vQhMjKSJUuW1HpE5de5SsnJyfz1r38lLS2N1q1b88QTT5xx7Lhx45g1a1b15zNnzmTs2LFnPe+cOXNwd3dn7dq1/Oc//2HKlCnMmDGj+vknn3ySF198kaeffppt27bx4YcfEh4eXofvXkQckfv5DxERaRypqanEx8fz4Ycfcv/99/PSSy8xfvz4Op+jefPmfPzxx4SGhgIwbNgw3nrrrTOOvfXWW3nyySc5cOAAACtWrGDevHksXbr0jGOjo6N57bXXMJlMdOrUic2bN/Paa68xfvx4CgsL+c9//sMbb7zBmDFjAGjfvj0DBgyode5BgwaRlpZGcXExUVFRzJ8/n6SkpDp97yJifypKIuIwUlNT2bRpE/fffz9fffWVTUUhNTWVa6+9trokAezbt48OHTqccWyLFi24+uqrmT17Nlarlauvvvq0r/u9fv36YTKZqj9PSkri1VdfpbKyku3bt1NaWspll11W57y/WrRokc1fKyINR5feRMRhpKamMmLECEpKSjh+/LjN5+jXr98Zj9U0QXzcuHHMnj2bOXPmMG7cOJte08fHx6avExHHp6IkIg6hsLCQvXv3MnHiRN544w1uuukmtm7datM5/jgZ+1xF6YorrqCsrIzy8nKGDBlS47n/OBl89erVxMbG4ubmRmxsLD4+PixevLhOeUXE8enSm4g4hLS0NNzc3IiLi6Nnz55s2bKFoUOHsnbt2hovh9V0ju7du1c/duDAAY4dO1ZjUXJzc2P79u3Vf65Jeno6kyZNYsKECWzYsIHXX3+dV199FQBvb28ef/xxHnvsMTw9Pbnwwgs5fPgwW7du5c4776zl34CIOCIVJRFxCKmpqXTu3BkvLy8AXnnlFbZv386IESNYtGgRnp6etTpHp06d8Pb2rn5s48aN1XfS1SQwMPC857799ts5efIkffv2xc3NjYceeoi77767+vmnn34ad3d3nnnmGbKysmjVqtVZF7QUEedisv5+IRAREScycOBAEhIS+Pe//+0Sr3M+JpOJzz77zOm2ZxFxZpqjJCJO7c0338Tf35/NmzcbHaXB3HPPPfj7+xsdQ6RJ0oiSiDitzMzM6i1IWrduXavLc7YwekQpNzeXgoICAFq1aoWfn58hOUSaIhUlERERkRro0puIiIhIDVSURERERGqgoiQiIiJSAxUlERERkRqoKImIiIjUQEVJREREpAYqSiIiIiI1UFESERERqYGKkoiIiEgNVJREREREaqCiJCIiIlIDFSURERGRGqgoiYiIiNTg/wELKSg8R7zljgAAAABJRU5ErkJggg==\n"},"metadata":{}}]},{"metadata":{},"cell_type":"markdown","source":"Pretty weird, right? :)"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3 (ipykernel)","language":"python"},"language_info":{"name":"python","version":"3.7.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":5}