{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Weak lensing\n", "\n", "This example computes weak lensing maps (convergence and shear) for a redshift distribution of sources.\n", "The lensing is simulated by a line of sight integration of the matter fields." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup\n", "\n", "Simulate the matter fields, and use them to iteratively compute the lensing fields.\n", "\n", "To obtain the effective integrated lensing maps of a distribution of sources, the fields from each plane are collected and added with their respective weights.\n", "\n", "The setup for angular matter power spectra matches the definition from the [Matter shell definition](shells.ipynb) example.\n", "CAMB is also used further below to compute the theory lensing spectra." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import healpy as hp\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "# use the CAMB cosmology that generated the matter power spectra\n", "import camb\n", "from cosmology import Cosmology\n", "\n", "# almost all GLASS functionality is available from the `glass` namespace\n", "import glass\n", "import glass.ext.camb\n", "\n", "# creating a numpy random number generator for sampling\n", "rng = np.random.default_rng(seed=42)\n", "\n", "# cosmology for the simulation\n", "h = 0.7\n", "Oc = 0.25\n", "Ob = 0.05\n", "\n", "# basic parameters of the simulation\n", "nside = lmax = 256\n", "\n", "# set up CAMB parameters for matter angular power spectrum\n", "pars = camb.set_params(\n", " H0=100 * h,\n", " omch2=Oc * h**2,\n", " ombh2=Ob * h**2,\n", " NonLinear=camb.model.NonLinear_both,\n", ")\n", "\n", "# get the cosmology from CAMB\n", "cosmo = Cosmology.from_camb(pars)\n", "\n", "# shells of 200 Mpc in comoving distance spacing\n", "zb = glass.distance_grid(cosmo, 0.0, 1.0, dx=200.0)\n", "\n", "# linear radial window functions\n", "shells = glass.linear_windows(zb)\n", "\n", "# compute the angular matter power spectra of the shells with CAMB\n", "cls = glass.ext.camb.matter_cls(pars, lmax, shells)\n", "\n", "# apply discretisation to the full set of spectra:\n", "# - HEALPix pixel window function (`nside=nside`)\n", "# - maximum angular mode number (`lmax=lmax`)\n", "# - number of correlated shells (`ncorr=3`)\n", "cls = glass.discretized_cls(cls, nside=nside, lmax=lmax, ncorr=3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Matter" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# set up lognormal fields for simulation\n", "fields = glass.lognormal_fields(shells)\n", "\n", "# compute Gaussian spectra for lognormal fields from discretised spectra\n", "gls = glass.solve_gaussian_spectra(fields, cls)\n", "\n", "# generator for lognormal matter fields\n", "matter = glass.generate(fields, gls, nside, ncorr=3, rng=rng)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Lensing" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# this will compute the convergence field iteratively\n", "convergence = glass.MultiPlaneConvergence(cosmo)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Galaxy" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# localised redshift distribution\n", "# the actual density per arcmin2 does not matter here, it is never used\n", "z = np.linspace(0.0, 1.0, 101)\n", "dndz = np.exp(-((z - 0.5) ** 2) / (0.1) ** 2)\n", "\n", "# distribute dN/dz over the radial window functions\n", "ngal = glass.partition(z, dndz, shells)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulation\n", "\n", "The simulation is then straightforward: compute the lensing fields for each shell and add them up using the galaxy densities as weights." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# the integrated convergence and shear field over the redshift distribution\n", "kappa_bar = np.zeros(12 * nside**2)\n", "gamm1_bar = np.zeros(12 * nside**2)\n", "gamm2_bar = np.zeros(12 * nside**2)\n", "\n", "# main loop to simulate the matter fields iterative\n", "for i, delta_i in enumerate(matter):\n", " # add lensing plane from the window function of this shell\n", " convergence.add_window(delta_i, shells[i])\n", "\n", " # get convergence field\n", " kappa_i = convergence.kappa\n", "\n", " # compute shear field\n", " gamm1_i, gamm2_i = glass.shear_from_convergence(kappa_i)\n", "\n", " # add to mean fields using the galaxy number density as weight\n", " kappa_bar += ngal[i] * kappa_i\n", " gamm1_bar += ngal[i] * gamm1_i\n", " gamm2_bar += ngal[i] * gamm2_i\n", "\n", "# normalise mean fields by the total galaxy number density\n", "kappa_bar /= ngal.sum()\n", "gamm1_bar /= ngal.sum()\n", "gamm2_bar /= ngal.sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Analysis\n", "\n", "To make sure the simulation works, compute the angular power spectrum of the simulated convergence field, and compare with the expectation (from CAMB) for the given redshift distribution of sources.\n", "\n", "We are not doing the modelling very carefully here, so a bit of discrepancy is to be expected." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAG5CAYAAAB1OMuOAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAe1xJREFUeJzt3Qd4k+XXBvC7ZQ9l7yEKyAYZguwtojI/EVCGiqjIHv4FF4oKisoSFyCKuEBQUUGmIIggIFORKcreUPZq813307zxbUhLG5Jm3b/ris1q+gbb5OQ85zknyuFwOCAiIiIiLtH/nRURERERUoAkIiIi4kYBkoiIiIgbBUgiIiIibhQgiYiIiLhRgCQiIiLiRgGSiIiIiBsFSCIiIiJu0rpfIckTFxeH/fv344YbbkBUVFSgD0dERESSgf2xT58+jYIFCyI6OvE8kQIkLzE4KlKkSKAPQ0RERLywZ88eFC5cONHbFSB5iZkj6x/4xhtvDPThiIiISDKcOnXKJDis9/HEKEBKoXfeececYmNjzWUGRwqQREREQsu1ymOiNKzW+wg0W7ZsiImJUYAkIiISZu/f2sUmIiIi4kYBkoiIiIgbBUgiIiIibhQgiYiIiLhRgCQiIiLiRgGSiIiIiBsFSCIiIiJuFCCJiIiIuFGAJCIiEoQeeughtG7d2u8/58UXX8Rtt90WNI8TLBQgiYiIBKGxY8fi448/RrCO6fj2228TXDdo0CAsWrQI4UKz2CTZOH8uTZo0gT4MEZGIwHEYoSRr1qzmFC6UQZJriouLQ//+/ZEpUyY8/fTTgT4cEZGwMmPGDFSoUMG8xubKlQtNmjTB2bNnr1pia9CgAXr37o1+/fohR44cyJcvHyZOnGju+/DDD5vp9CVKlMCPP/7o+h5moLJnz57g5zHzk9Sg1tWrV6Np06bInTu3CdLq16+PtWvXum4vVqyY+dqmTRvzONZl9yU2vncMGzYMhQsXRoYMGcxtc+fOdd3+zz//mO//+uuv0bBhQ2TOnBmVKlXCihUrEAwUIEmSOMuYwdGYMWNw+fJlvPHGG9iyZUugD0tEJCwcOHAAHTt2xCOPPIK//voLS5YsQdu2bc1rrydTpkwxgcuqVatMsNSjRw+0a9cOtWrVMkHMnXfeic6dO+PcuXNeH9Pp06fRtWtX/PLLL1i5ciVKliyJu+++21xvBVD00UcfmeO3LntaInzrrbfw5ptvYuPGjWjWrBlatmyJ7du3J7jfs88+a5bn1q9fj1tvvdX8e1y5cgWBpiW2FHrnnXfMictNkeC5557DuHHjXJf5Rzty5EhMnjw5oMclInIt1apVw8GDB1P95+bPnx9r1qxJ1n0ZYDAYYFB00003meuYTUoMMyx8XaYhQ4bgtddeMwFT9+7dzXUvvPAC3nvvPROQ3HHHHV4df6NGjRJcnjBhgslC/fzzz7j33nuRJ08ecz2v43NNDAMjrjp06NDBXH799dexePFi84Gb76MWBkf33HOPOf/SSy+hXLly2LFjB0qXLo1AUoCUQj179jSnU6dOhdz6cEoNHz7cnCxM/54/fx5Tp041v8RFihQJ6PGJiCSFwdG+ffsQzBjwNG7c2ARFzLAwA3TfffeZJTRPKlas6DrPmlAuydkDKi670eHDh70+pkOHDpkgjNksPg4TAsxI7d69O9mPwffI/fv3o3bt2gmu5+UNGzYk+pwKFCjgOn4FSBKUmDVi2tMyfvx480nn1VdfNZ92mDblpwARkWCVVHYjWH4ug5wFCxbg119/xfz58/H222+b197ffvvN4/3TpUuX4DJreOzXWbVFrP+h6Ojoq5brWC6RFC6vHTt2zCyRMavF+qGaNWvi0qVL8Iekjj+QFCDJVT788EP07dvXdZlpUWbNjhw5glGjRpksEgsD+QmDqV0RkWCU3GWuQGNQwMwKT1wiY1DyzTff+OSxuRzG2iEWcmfJksVcx1qfpCxfvhzvvvuuqTuiPXv24OjRo1cFNUmVmtx4440oWLCgeSwWedsfu3r16ggFKtKWBL744gvXWjbxj/V///uf6w/t0UcfNeeZbrXXJomISMoxU8RSBgZzXMLiji5+GC1TpoxPHr9GjRpmd9gzzzyDnTt34vPPP79mbyUWZU+dOtUUjfP4HnzwQVNiYceda+x5xGXMEydOeHycp556ynzAnjZtGrZu3YrBgweb4Mz+ATyYKUASl1mzZpndD1Y6dsCAAWbbph2L6dKmjU88MhVs7WoQEZGUY6Zl6dKlJlvDHVzMzLOEoXnz5j55/Jw5c+LTTz/FnDlzTK0SPwS7v657WkU4ceIEqlSpYt4T+vTpg7x58ya4D4+RS4OsRa1cubLHx+H38X1k4MCB5mdzi/93331nArBQEOVIbC+hJMkq0o6JiTG/4KGOv+jcnWCtMT/++ONmJ4SnXhnst2F9AuG2fwZNIiIi4fT+rQySYNmyZWjVqpUrOOrUqZNZf06skRi3bVq3sSbp4sWLqXq8IiIi/qYAKcJx3Zv9J1h4bXVGZfMv7nxIDLde8n7EnW1sXCYiIhJOFCBFsE2bNpm+G1Yd0V133WXWp60ao6SwQZmFRXjB0PVURETEVxQgRaht27aZWTvHjx83l7kNc+bMmabfRXI71HJeEP3999/46quv/Hq8IiIiqUkBUgRiZ1kGN+yWam0D/f77781W0JSwZ5HY7l71/iIiEi4UIEUgzlJj4y+rzT0nP3MKdEpx+rLV8Itzf7iNVEREJBwoQIpAv//+u+s8g5rEZv5cC3ey2bNII0aM8MnxiYiIBJoCpAjDZTB2R6XChQubVvDXo2XLlihbtqyrhTxbBoiIiIQ6BUgRhi3srcJsX7SyZzsA9kWyKIskIiLhQAFShNm8ebPrvK9m/XTs2NEMVyTWM11rEKKIiESuhx56CK1bt0awU4AUYazlNV8GSJzqbB83wh1tIiISHpYsWWJqTk+ePJmi7/vnn3/M97l/aB47duw1B+YGAwVIEcYfARI98sgjyJMnjznPnkg7duxAMNdhsXdTXFxcoA9FRCTiZMuWDdmzZ0ewU4AUYfwVILGHUv/+/c15Bh5sJRCs+vXrh+LFi6Ndu3aBPhQREfOayfrNm2++GZkyZTLtV2bMmGE+zLFnHSceWH3mWEPKDTYvvPBCguzO7NmzUbFiRWTMmBF33HEH/vjjjwQ/45dffkHdunXN4xcpUgR9+vTB2bNnXbdzpibrSXkbGwaXKFECH374ockCsaULccczfxaXyGju3LmoU6eOCXZy5cplBp7v3LnT9Zh8PlS5cmXzfQ0aNPC4xMafzePJmzevOX4+5urVq123W89x0aJFpkkx329q1aqFrVu3wq8ckiLjx493lClTxnHrrbfyt9URExPjCCWFChUyx50zZ05HXFycTx/75MmTjhtvvNE8fvr06R379u1zBJvY2FhHtmzZzDHytHv37kAfkohEuFdeecVRunRpx9y5cx07d+50fPTRR44MGTI4lixZ4ti7d68jR44cjjFjxpj7tmvXzlG9enXH5cuXzeXFixeb1zK+L82fP9+xceNGx7333usoVqyY49KlS+Y+O3bscGTJksUxevRox7Zt2xzLly93VK5c2fHQQw+5juH+++93FClSxPH111+bY1i4cKHjyy+/dFy5csUxc+ZM8zO2bt3qOHDggHmtpxkzZpjbtm/f7li3bp2jRYsWjgoVKpjXWVq1apX5Pj4Wv+/YsWPm+q5duzpatWrl+tl9+vRxFCxY0DFnzhzHn3/+aW7nc7bubz3HGjVqmH8T3qdu3bqOWrVqefXvzfft5Lx/K0DyUnL/gYPxmHmqXbu2X37G4MGDXT9j4MCBjmDDP3Dr+Hh6++23A31IIuIvVavyU2Hqn/hzk+nChQuOzJkzO3799dcE13fr1s3RsWNHc3769OmOjBkzmtdXBjoMcixW8MBgxsLAIlOmTI5p06a5Huuxxx5L8PjLli1zREdHO86fP+96XVywYIHHY1zs/BknTpxI8rkcOXLE3G/Tpk3m8q5du8xlBk929gDpzJkzjnTp0jk+++wz1+0M7BgwjRw5MsHPZ6BlmT17trmOx++v9+9rTyWVsLFlyxa/LK+5L1+NGTMGFy5cwPvvv49nnnkGOXPmRLBYtWpVgsvffPMNevXqFbDjERE/OniQs5UQzFivee7cOTMb0+7SpUtmaYpYDsDXKm6Aee+991CyZMmrHqdmzZqu83zNLVWqlKukYsOGDWbawWeffea6DxMkXNrbtWuXGVyeJk0aM5MzJbZv326W+n777TccPXrUVde5e/dulC9fPlmPwSW5y5cvo3bt2gk2/nBKg70khLiEaClQoID5evjwYRQtWhT+oAApgvir/sguX758pmD73XffNevb48ePd62VBwP7ujb9/PPPOHbsmFk/F5Ewkz9/0P/cM2fOmK+sISpUqFCC26zh4QygOAGBQQyDkpTiz3j88cdNnY87Bhfebqpp0aKFafEyceJE03SYARIDIwZ3/sDAycKaJPLnZhsFSBEkNQIkeuqpp/DBBx8gNjYW48aNw8CBA5ElSxYEYwaJx8gXpi5dugTsmETET9asQbDjJAIGQsy6JJbB4Wsom/Kyz9zdd9+Ne+65B40aNUpwn5UrV7oyKSdOnMC2bdtcr/NVqlQxPfBYeO1JhQoVTKDBD4wsCneXPn161+ulhR8sWSTN4IjF31Yh+LW+zx03zPB+nMRg9dNjRokfZrkiEUjaxRZBUitAKlasmGkeaf0R8Q8oGPCPbt26dQn+cImpaxGRQOCgcPaR4y7gKVOmmCWntWvX4u233zaX+QFu8uTJZnmMy3D8ANq1a1cTBNkNGzbM7PLi7jXuEsudO7drpxh3p/3666+mnIA9iZiFmjVrlqu8gK/ZfExm/7/99luz7LZkyRJMnz7d3M7AhRmbH374wUxjYEaKO9qYeZ8wYYLJQP30008YMGBAgmPirjTumuNut0OHDiEmJuaq588Pzz169DDPi/djINe9e3eTNevWrRsCKsXVTRKyRdolS5Y0x8yCQGuXgb/88ccfrkLo4sWL+/3nJcfvv//uOqYOHTo48ubNa86zmPHs2bOBPjwRiVDcUcxdaqVKlTIFy3ny5HE0a9bM7NjKly+fY/jw4QkKmKtWrWp2ndkLmL///ntHuXLlzA5i7nLbsGFDgp/BHWVNmzZ1ZM2a1RR6V6xY0fHqq6+6bmexc//+/R0FChQwj1GiRAnH5MmTXbcPGzbMkT9/fkdUVJQpsiYWdXP3HHfc8fF4vDyWb775xvV9EydONLvjWBBev359j7vY+LN79+7tyJ07t3ksbiLi8SZVJM7Cb17HQnB/vX9H8T+BDdFC06lTp0yzK0bEN954I4Id+0ywdwTTqCz84ycUf+OnnYULF5rz/GTAXh6BxKJxflIhFpL/+eefruwWs0ih0PpeRMSOmR72KWJGKRSaL4bS+7eW2CIEU6pWMZs/l9fsnnzySdd57rwIpgJt7pCwB0RaZhMRETsFSBEiteqP3Hc4WLsyvv/+e1OEGAwF2mnTpsVtt92Gxo0bI2vWrK7ju3LlSkCPT0REgocCpAgRiACJgchjjz1mzjN7FchibRYVsvjP2rHBwkHuHOGOEGJ6eunSpQE7PhERb3B8BytltLzmewqQIkQgAiR69NFHTaBEDJD81R/jWlhzZS0x3n777a7r27Rp4zrP3RsiIiKkACnCAiQ2GkusF4Y/sHmYVevDbZ6BqvVxrz+yMINkNR9jgKQ9CyIiQgqQIgCbdFlTjxkc2XsARUqxtr1BpD2DxB0MrEWiPXv2pMruPhERCX4KkCLAv//+a2ajpfbymn2NvHTp0uY8O7Vye32gAiS2OmDnWjv7Mpt2s4mICClAigCBqj+ysAOr1X/I6keUmtj59Z9//jHnq1at6qqJsrRs2dI110cBkoiIkAKkCBDoAIk464zZG2L7fGtAY2rXH9mX1yz58+d3TcLmTjfOMBIRkcimACkCBEOAxC2oDzzwgDl/+vRpfP755wEv0LazN43UbjYREVGAFGEBklULlCLc2cUhg9u3Azt2AFyu2reP29KA48cZ8QDnzwNstJjELjD7Mts777yTajvGEivQjuQAiUN7N27cGOjDEBEJWgmLMSTsMAixAqQiRYq4Oke7sC/RgQPA/v3xQY91sl/m+bNnk/9DWePDrfP8mikTkC8fUKAAqhQogI8LFMDvBw7gwMaN+OO991DhzjvNbciSBf56/lYGiZOnb775Zo/3K1myJMqVK2cKyFesWIEDBw6gAI8rDLEhptVcjoESu4qLiEhCCpDCHHsPnTx58r/ltcOHgU8/BWbMAHbujL/sa8wkWWM7mF3iz9i0yVzs6jwZPXv+9z033BAfKPFUpAhw661AqVLxp5Iluf3M6x18LNK2skdWMbYn3M1m7bD77rvv8PjjjyMcffDBB67s3YIFCxQgiYh4oAApzDF7xP/JHKjxCpfHOBstJTPHGLjwewoWZDUzt6QBly/HP4b9a2LXMfN08GD8+aQwkOIpsQJpBk1WwGSdGEQVLQpERydreS2x+iN7gPTKK6+4drOFY4B0/vx5E/xZrP5YIiKSkAKkcPbHH8jxyivYCyAfL//9d8LbCxeOPzH4YRBkBUL28wyQrhezFaxV4lLewYP4YtQorPvxR3AB666KFVEmR47423hikOTJnj3xp4ULE16fMWN8hqlcOYCZEOvEZb1kFmhbKleujKJFi5qhuj/99BNiYmKQLVs2hJN58+Yl2EGoAElExDMFSOHmxAngiy+Ajz4C1qzBVYsnDHq6dAEeeig+C5MamHXKlSv+VL48bi9WDA8wqAHw9unT2L52rRmBYvDNm0XgfON2PzmXChNgA0wu3/H05Zf/Xc9sV6VKuG3zZrQHsJ5LbFWqXOMwo0yx9rhx43D58mXMmTMHHTt2RDBmBVeuXGl2BXLgbkp89dVXCS4rQBIR8SzKoeFTXjl16pTJLjDLwHEVARUbG59ZYVDEHVgXLya4mZdmAWj2xRfIdt998cXTAXbXXXeZbAbNnj3bzERLEn9Njx79L1jiUpx1nrVU11rCIxaMV6jwX5apWjWgYkXAFmQsXrwYjRo1MufbtWuH6dOnI5icPXvWjIs5ePAgOnXqhKlTp6ZoeS1v3rxX9aA6fvw4cjCLJyISAU4l8/1bAVKoB0gsgG7QgGmFq2+rUgXP7tyJ92NiEJUrF44ywAgSs2bNcm2tv/fee/H99997/2AMjhgorV+f8HTs2LW/l3PpGCxx+a16dVypUgX569bFsRMnzI4/Fnhn5DJekPjhhx/QokULcz46Oto0tSxevHiyvpftC+xjVSzctXfHHXf4/FhFREL5/Vt9kELd3LkJg6PcuYF+/YANGxDz008YHhOD4wFsEJmYe+65x7QdsDJI1igQr7ClQPnyQKdOwJtvxmfTjhzBtDffxD0AnuVSEjNHJUpc/b1sc8BC7vHjzdJj2vLlsfvMGcwHMOTMGfzBom22OQgSc/n/2ykuLg6jRo3yanmtadOmrvNaZhMRuZoCpBRig0MOO02s4WCqs9flPPdcfN+i0aPN0tGWLVtcNwVbgMR5aI899pg5zyTmhAkTfPsDoqKwZMcOzAEwHMCBcePiG12eOgUsXw68/XZ8LZaHxpmZL18Gw4dnAFR79dX4gnUWs7dtGx+A/fprfO1TgAMk+uijj5KVGbTvXuNy2pNPPum6TQGSiMjVFCClUM+ePc28LvvuqICyN3BkTQ2XjIJoxEhSHn30Udfg2EmTJuGiW+3U9bK2+LP4mkNqDe7Kq1UL6NWLQ+His28sbGfWafhwttRGHAu83THw5CDbp54CatcGuLuN89sGDgRmzozfgednO3bswE7WW7kFPu+++26Kdq9xma08M25OCpBERK6mACnU2Qtu3bpRB3uAxCGxbZmVAVfEjuDrr7/22WNfuHDBNUqDGb8bkmpXkD070LgxMGSICYKiDxzAo3fdBVbrjODGQAZX7uvUXJpbuRLgEhcL37k7kF26H3yQaUbO8khZv6lk+PHHH13nn3jiCdfOv/Hjx5tAKSn2YnMWnxcrVgzpuDSpAElExCMFSKHOnkFyGyMS7AES2Zd6kpMJSa7169fjijNA8WY5tP4DD+Bb5zLbK/Xrx2eZ/vgDmDgRePhhzy0SWEfFIbzMTrGlgBV4DR0an6Fy2z12PctrnGvHQMcKLj/55JNEv4/Bk1UEz+W1xo0bm8wdd8NZmalY7oQUEREXBUjhFCAlkkHKnDmzqyA62NSrV89keOiXX37BJudIkuuVkgaRiRWRWxkadtV2sJcTm1E++igweTLA+i6OMGHgwcwTgyi2EXD/f/PTT8CwYayKjg+YeCwDBsQv1zlHoCQ3I8YWBMQZcRUqVMCgQYNct7/11luJBjkMrOzLa1bmqJQzyOPSJkeyiIjIfxQghekSG99Q/3Z2zi5durTZEh6MWB/EbIjlvffe88nj2keMeJNBypkzpxnoSrt27fIcuHHH4L33xtcuLVkCxMQwMgPGjgXat48v7LZjAMPbWUTPpcW8eeOLxLt3B5gB4v+vRLpuLFu2zLWMxh5SVl1Vw4YNzXXbt29PMEIksd1rVtbJHiAR2wWIiMh/gvNdU647g8Q3TG4DD+blNUvnzp2RxXnsbHx4OrFxI15kkNKnT4+KbAbpBatPk5VFuiZmZlgo36dPfFdvjkZhZobDgTnXjRkod6z/mTQJ6NoVYD8jBlU8z2DHtlPOvrzGAMnyFIvGnd7kDrtr7F7j8pqnAEl1SCIiCSlACtMapFCoP7KwYdeDLG42CbEz+JQBxXU4efKk6w2fk+oZJHmjVatWCZoseoXDdPnc3n8/voaJW/JnzQK4PMbmjO5dzdlzidkk/uw8eQCOOpk5E0tmzzY3MxPYpEmTBMFSOWfg9euvv5qTHQMrdt92X16jWzns10kBkohIQgqQwilAstXAhFKARPZlNhZrX0+D9zVr1lxX/ZGFdVvVmBFyFn1zqe26cR5dy5bAG2+whXX8shxri6w6JXsdGZdPmYm67z4s27oVXCh7rnhx5LQFVVxqs9civcHHTWT32v3335/gNmWQREQSpwApXGqQMmdmeiFkAyRmemqyrxCYaPkDS5cuDViBtp19NIfXWaSk8P8ba52efx6YPz9+PMoPP8QPE7bNR8sM4D4AL7HZJWuXGGRNm8Y1NDO0tiDbDDhHuFj1RPbda6ypsmbMWXLnzm2uJwVIIiIJKUAKlwxSIlv87du5g519y/8LL7zgdRbpegu0UzVAcsfBuffcEz94+NAhdnjE/GLFcNh+HzbUZODToQObSSH9k09i5L33IsrZldwaP8K+SYktr7lnkfbt23fVEFsRkUimAClcAiTb0gy3e1sZAQZHnt4Yg1H79u1RsmRJc54ZJPexGinNIHEIob3OxhvcAWg9BtsQsOeQJyyIZzBy+PBhsxTHLNhvv/2G3bt3e//D06XDpQYNcN+xY2B+qNWNN8LBILJAgf/uw9EpH36IBydMwK6oKAzjcU6ebI4jsd1rdtrJJiLimQKkUGd96rcFSBz8ao3tCIXlNQsDuVc5+8xpyJAhrp14ybV//36TDSHWD11vewPW+FhZJB4LC6Rr1KhhRnXcfPPNyJs3r9mBx55JWbNmRb58+XDLLbeYPkV33HEHihcvbobxemvFihVmVx87HGVllohduvfuje+vxIaVtszhTQ4HnucS5eXLuFCpEvJ//TWyJbK8ZlEdkoiIZwqQQhmDB2vERIjuYHP3f//3f665aRs2bMAXX3wRsPojT9v9Ob6ES3h//vmnCUSZUTp37lyi38tu3o888kiimadrsWfRmjdvHn+GQR/7H7FhJZfhPvuM29ngsAWDRQ8exOhLl7AfwHe5ciEdC9c9LFkqQBIR8UwBUiizvzHbMkihHCAx4/Paa6+5Lj///PO4xLlnAag/sgda97IhpA2zRsweMYvEbBKzSszStGjRAh06dDCDeK2fz+Wuxx57zKuaKnuAdOedd3ou8n7gARYcIWrvXnxWpQrW2W8GUJuF3RzQW6ECMG5c/NgUJwVIIiKJcIhXYmJi+G5nvgbMgQN8y40/tWzpuvrhhx82x8bTmjVrHKGoSZMmrufw9ttve/V9e/bs8ekxnThxwnH27FlHXFxcsu5/8OBBR548eVzH89FHH6Xo5+3fv9/1vVWrVk3W96xfv97cvwLgGAs4TkZF/fc7Yp0yZnQ4OnVyOJYtc1w4f94RHR1tvqdy5copOj4RkXB+/1YGKQy7aNszSCwyDkUjRoxwnX/55ZeTtcOKNUJWDyTOKytUqJBPjyl79uxmrh3rkpKD9UgTJkxwXe7Tp49ZlkuuefPmeeyenZRKlSqhadOm4GCUvqzj6tIFmDIFqF37vzuxQzebcdatiwy1amFgrlzI4CzSvlaWi0uMbORpjT0REQlXCpDCrIs23+CsAKlo0aKuER6hhgXW1s4rLlGN5vyya+BUenbRJi5vJTeQ8SfWLz3EnkaAKbbu2rVrokNlk1V/lAxDhw5FhgwZzKlb794Ag6RffgH+/BPo1y9BfyWsW4eRR45gD4DBZ8/igK3JprujR4+ibt26ZjQMg7AYNrkUEQlTCpDCLIN08OBB1xtXqNUfuXvllVfM7jCrQ/S1Cp39UaDtC2PHjsVNN93kal+QnGCPQdR8No50jmJhjVNy1a5d27QZYCG5VfBulC0bPyjXGmdiq9HKwy7dAPKzWSfHm7DLt1s2ac6cOTjFtgIAli9fbmquGDSJiIQjBUihzL7s5AyQQrlA2x37D7HY2cq+DB8+PNULtH2B/ZimTJniymg9++yz2LSJi2BJB3snnMXUzNaw4WdKsP8VWwx4lDEjJwQDv/1mAqHt1arhsvOmaGa3ON6ERd1cluOQXmerBQZIdmvXrkX9+vVx4MCBFB2biEgoUIAUZkts4RQgWR21MzlnzHFG27///pusDJI1Qy1YMJAYOHCgOc9deVymsnpVXWt5Lbn1RynGgO2OO7B35Egwv/UyA1HbPD+TRWrblr9IiH3vPSxxHhMDPmu0yebNm82yW1L/X0REQpECpDBbYgu3AIlvxH379nUFFqyv8eTy5ctYt26dK3tizRgLJiw2Z0sAq8fTiy++mOh9OSbE0qxZM78eF7f6Mwf0AoAH69SJH3PiPE5j2zakefJJrI+JwRDOhGvcGMuWLUOxYsXMzTt37kSdOnXUiVtEwooCpFAW5ktslqeffho5nIXFn3zyiamvccfrLnB3VpDVH9llzJjR7ACzRr+MHDnSjC9xx7oeKxvGjtyFCxf263Fxxx+7gNMfO3bED8rduNH0VjINKZ3yA+Ai5/tz5uCW997D8pkzXWNY9u7di3r16l1z6VBEJOwCJM6VktBZYuOkdp7CAbfXDx482LVL75lnnkmy/ihYAyRrG/6wYcNcbQm6dOli6qvsFixY4Npu77flNRvWRlkNI9mGwASaXH7jz+ZIk9WrMS9bNjPuhNJxafDNN1Gwdm2sq1cPTZytJA4dOmSaWVqF3CIiEREgsbDUnkJftGiRmUslwbPExt1rVsFsuGSPLL1793b1Nfr+++/NLqpQKND25KmnnjI7zYiDbQcMGJD69UdurACJgRnbJdjtzZ8fd8XEgLmiGXnzAhkyuPopZZ40CfP//hsz8+RBMecuyuTs0hMRCZsAiUsDjz/+uPnKIZysBWHjOwmeJbZwXF6zsFDbXrPDjJK9qaG1JMW2AJUrV0Yw4zFyqdBa1po0aZIJ+qyskhUgsYeVFUj5W1IjR6x6qL8B/Pnkk4zqAAZ1HHPCDNSlS2h75Ai287kAmD5ypNez50REQi5A4qfK119/3dSDcD4Waydatmzp36OTFGWQwjlAIjZctN7I+fs3e/Zsc/7s2bOm549Vs2Ptegtmt9xyS4JMC9sZMKhYv369aYxJ7DPEZo+BDpCsf2e6++67WbQEvPUW1+OAIUOAG24wt7ERQTf2njx3DtuYXd67N1WOXUQkoAHS22+/baaScxmHNSD/+9//TG8XCZ4apHAPkNgL6NVXX3VdHjJkiGmoyH48zLwEe/2Ru27dupnhtsSgiBla++61lHTP9leAxFYECxcuNOc5nDdB48k8eQD2pmKg9OKLiLvxRnN1ejar3LgRjhIlAO5APHgw1Z6HiEiqB0jTpk0zO4XYvI7LbFzG0LbeAIuwDBK1bdvWFQTx9/Hzzz8Pqfoj9+LoiRMnuorpv/nmG9MxPLW299uVLFnSY4DE7fzM0FkBW3S0h5cMtlQYOhTR//6LhTVrwio5j2Ix97hxTJdxKyLgHAMjIhKWu9i4RZnLGB07dkzwaV6CpwaJtStFihRBOGJQ8dprr7kuW8u9llDKIHkaaGuNieH2eS7DpRb77wwDJKu+y9492yyvJSV7dlT54QdUuuEGvA7gnHU9B9uOHBkfKI0axbSU356HiEjAd7H99NNP2sUWRBmkC2nTmh1RVLp06aAY1OovDRs2NNvJiR2cv/32W3M+c+bMKMt5YyGmTZs2roG2qb17zdMyGwf+WkXWVv0RC8utf/OksEFnt6efBpsy3MzvZ2bKqqPi6BR2E2dbgM8+c40wEREJq11sHAGhXWxBEiClSYNt//zjqsMJ1+U1O3sWyVKlSpUUzywLxoG2gQqQrKaPxA9D3JhhfSjibjr2o0oOdj5nvRJLze/dvh2bZsyIn/1mBe2sWerUifNg2C/EP09GROQ6aRdbOCyxcXltyxbX1ZEQILEGrkOHDgmuC7XlNTvON5sxY4bJgLVv3z5Z2Rp/F2rbC8avubxmw/YFfI2wDHr7bbZA53RbwP68OBqGWWgGgxs2mHqnmjVrmuL1K1eu+OIpiYgEdhdbv3798NFHH+H3339PcgCn+CmDFAE72BKbbWbPGIVSgbYnHLDLdgVffvmlWdIKdICUovojN4899phrVtv8+fOxZMkS4LbbgHnzeEX8ecu8eXBUrox/6tfHvpUrMXnyZEyfPt0XT0lEJLC72NivZd++fSbDxGUODuTkp2AWcVsN8MSPAVKE7GBzx6G0PXr0cBUZ169fP9CHFNLsARIH/y5evNicZ/G2NWQ3udKnT4+XXnrJdblnz56m+74p/m7aFPj9d2DqVDiKFjW3Rzkc6OxwgK8oXDx97/XXEzQCFRFJbVEOP7wKcZYTg6mNGzea4ZXhOHqA86ayZctmMmpcHgmI9Ok5xt58Gq8YG2v+rZlROXfunGsgarhjHyR2pWZQyNo48R5r2BhoWkN/Law9fP/99736f8PZc1YTTypXrpwZG9OpUydcunQJne67D6V/+gnPAYgfRxyPA3NOPP00yrLPkqfWAiIifn7/TlGAtHnzZnzxxRcYOHBgsgs2w1XAAyQGRgyQOD+rdm1kWrPGLG8yUOD/JxFvVKxY0QTadrNmzfK63pCZqFatWmHPnj0JrufrB/9+uAuR8qZNiwWNG6PcokVIY68/YtDLXkp+Xj7ly+Dw4cOxf/9+kwm3xsCISOS+f6foo9mIESNMZshTcMRPnVtshcKSelv8z0VHu2q/ImV5Tfy/zGYtlTVu3Njrx+NS/M6dO80SvX2uHFsJWMERG2XO+OknVJw7F44//sB8+6iYlSuBGjU4i4XtxuEvrJF67rnn8O677yboSyUikStFAdLKlSsT3dqfMWNGdO/e3QRRkroB0inbJ24FSOLLAKlBgwZm2e16cLn3/vvvNztf16xZg65du5rAy8pYcdBw3bp1zeW0pUrhrxEjwL1urqo6Jrk//JDtvgEu1zN76mO//vqr6zw3m4iIpChA2rt3rymMTcwTTzyB7777zhfHJSnoon3MtnNQAZL4MkBK6e61a+Est48//tgsuXHGG8fEWLvdLNwtuzp7dlRki4DoaMQ5h+Hi1ClgwACgUiVgwQKfHhfn+Vm0RC0iKQ6Q2CX3wAGWT3rGPjTslySpm0E6bDuvAEmCOUCysJEkl+4yWF22bW644QbzYYt50bfi4vDaww9zsu9/jSa5Y5P9lNq0AXbv9nmAxFIBFpiLSGRLUYBUr1498+kv0QeLjr5qB0ww4+467qphcz4uHYbUtmJbULSfn6wTeYMTSQn77w8H2NqH2KYm7nSzdmK+8cknODNmDAdCxhdtWzhihqNl3nrrupbdjh8/jn/Y3duJr2H2yyISmVIUIA0aNMhMH0+siHHFihWpOmDzenDW1Pjx4029AXft8CtrrEJxiW3PsWPmK0dVXG+9iEQ27uzgqJAcOXIk6GOU2goWLIgHH3zQVdDN5pFmJ9vy5QAb1ObL998HhUGD4seWrFjh1c9av379VdfZWxOISGSKTmn9AHd5PPnkk2jatKkZErp7927zCYxbgTmG5IEHHkCo4DgDflq8fPmyOTHtH4oZpGOXLpmvWl4TXxgzZgyOHTuGjh07BvQ42E7Enu0140fYE6lLF66Dsfvkf8tuGzcCtWqxhTdTQl4vr1lUhyQiKe7Axp1q3BLL/gFt27bFzTffjDx58piJ5BUqVED//v19cmBLly5FixYtzCdJTqa3JrbbvfPOO6bAkzvoatSoYQo+k4vHzIxY0aJFzc9o0qQJihcvjlAMkKxzCpDEV/g3F2js3m0N7eWSFztxu7DVyPjx8ctulSv/d/3EiUDp0jg1fjya3Xmnmet2rWV/BUgi4olXLWrr1KljghG+iEydOhUffPCB2SY7d+5c1/bd63X27FnThZdBkCfsqzJgwAAMHTrUvMDxvs2aNcNhW6+U2267zbzIup/YDI4jU3744QfzwssxKTx+BmWhuMRmBUiBqhcR8Rd23LY3nbwKl934wWjsWFZ3x1935Ahu7N0bgxcswK+TJ5uRR8kJkPjaxTpKUoAkIv9N+vRC6dKlzckfmjdvbk6JGTVqlMlmPcwdLoAZhTB79mxTqzB48OBEawssX331lWlZwJ15dM8995gaJBaie8JGjPZBvOzEGSwZJCtUsp6LSLiwz4BLtC6IA4vZn+3//g9gBvurr8zVDfkawAa3w4dj8//9H8raB+Q6nT592jVTkh+y+MGJO3E525CjV6yASUQiT0j+9XOGE4uquSxm4QsZL7NQPDk4gJNZI6bfuaWXy4ZJ7QBjA0wWsFonfn+wLbFxe7RIOOHfpBWksIt/kgoVAqZPRxwzw84lQjYReDEuDhlr1UKcrRmkZcOGDa7dqxy0zR2txHmGrK8UkcgVkgHS0aNHTVCTz9rJ4sTLBw8eTNZjcLApe7xwFAK7+bL+KKl5U0OGDDF1V9bJfbZUMARIARuaK+InrC+0mtMmtz/RzltvRRnOVuNGDOd1t5w/j6g6dYC+fRMsT9vrj+wBEmmZTSSyhWSA5CusTWAqnan7cePGJVmYyoZ2DEDsp2CrQVIGScJ5mY3Z3r///vua92d2mWXZzwLoWaMG1jivj2KmiINvy5UD5s411ylAEpGwCpA43DJNmjQ4dOhQgut5OX/+/IgIHmqQAh60ifgBm7mmpD8R571ZWj7/PN7r0gWDuGxmXcmlM9Y3duqEv507X/l6wkBMAZKIXHeRNj/Nbdy40ewaYzGjXVJLVb7A3SbsycRtv61btzbX8Rh4uVevXogIqkGSCA2QrL/5xNiHzfJ1gsvppefMwTdHj4ItbhtbN372GWYC4CvGlrJlzXKefdOJAiSRyOZVgMTt/F26dDG1QO64TOWLOUZnzpxJMNdt165dZlcad2qxdxG3+HMqeLVq1cwMODa3Y2sAa1db2NMSm0RggHStQm1+ULKWzQoVKuTKKLPRZOfOncFtHb0yZ8aYtGmR5tQp5GHLEACruSv18GFkyZvX9Hbj6w0DJBZwW69pLNpm37Vg6BElIqnA4YUSJUo4nnzyScfBgwcd/rJ48WJuLbnq1LVrV9d93n77bUfRokUd6dOnd1SvXt2xcuVKR2qJiYkxx8OvAVGrFisqzCkN4MiQIUNgjkPEzy5evOhImzat+XurUKFCkvfdtm2b67WiZcuWruvj4uIc//d//+e67e4qVRzbK1Vy/Q2ZU+7cDse0aY577rnHdb89e/aY723fvr25/MQTTyT58/fu3et49tlnHb/88ovPnr+IBOb926sA6YYbbnDs2LHDEYnGjx/vKFOmjOPWW28NbIDkfHG/EBVljiM3X9xFwlTZsmXN7zk/DF2+fDnR+33++eeu4Oall15KcNuJEyccxYoVc91+4403OtoBjiP2IAlwbLj1Vkdu533mzZvnmD17tut7cuTIkeRxNm/e3NwvZ86cjgsXLvjs+YtI6gdIXhVp33fffaZvUCTq2bOnSb2vXr06KJbYtMVfImmZjT3Q7EvvSdUfcfndLnv27KYDf7p06VzNXtlSko98pVUr1/0qbtsGVh/d5+yTxOV8CxtJss2HJwcOHMC8efPMec6n5A5ZEYmwGqTx48ejXbt2WLZsmZm/Zr3gWPqwq62kSpH2aWeTO9UfSbgHSOx+bxVqJ9bB371A2x3rFV9//fUEQU/O0qWR9ptvOL8I4CaPY8dMbRJ/2jfPPecaBm1hfRLHGLmbziaVtg0rDK483U9EwjhA+uKLLzB//nyz64OZJHvRIs8rQEq9AEkZJInEQu3/41iRaxRouzeStfTr1w+LFy/G999/by6zWSz4GtahA9CwIS537450ztvaXLqE2gB6APj6GgESXxftGCCJSOjyaont2WefxUsvvWRSzRz2yhcM65ScRm5ynZg1cguQlEGSSO+FxKU3a0aip+yR/UPcxx9/bLJJWbNmRe/evf+7MV8+pJs1C71y5sQx51V5AdMOgOFPLsDjaxyv++233xJcl9QsSPf+bVeuWD2/RSSkAyTWAbRv316DHAPlwgV+XDZnlUGSSMBxI+x/llSAlFT9kTu2C+HcxpMnT6JmzZoJb4yKwo7bbze1Sd/Yru7A7BWATIsXX/V4X375ZZJz3hLD7ytQoIA5XgVJIsHFqwiH/YdY7CjB00VbGSQJZ6xztIZJb9u2zXxIS2n9kTt+wGMHbU/YVZt9+tsySLrvPsRlz26uZ1elJ2fPBp54IkEvMvvymtWNm4Xa+/btS/IYWFfFIIrB1DWH8YpI8NcgsWnayJEjzY4NDnp1L9IeNWqUr45PPNGgWonQZbZNmzaZTMv27dsTLLt5EyAl5YknnjB1SmwMedcnnyDqxAn8WKQImltF2B98ACxaBEydij+yZnUFN8xG1atXz9WFm8tshQsXTvTnMIjCNWqbRCSEAiS+SJnCxmR0tg0377zzjjn5olu419RFWyKQex2S/TILtK0AKakC7ZQs6dkDLmTKhAElS+LbrVvBj39ZeB3bDdSujSM1a5oXUi6QdezYEXnycA9cPGaG7r333kR/DtsGWFS/KRIGARI/WUUq9kHiicWg2bJlC8xBaFCtRCD3nWz3339/ggLt06dPJ6v+yFs333ILJmzdip8YoFWtivQMoOLi0HD5cqwE0CUqyhyTPSt0rUJt9wySiIR4gDRs2LAkd4g8//zz13NMci0aVCsRKKmdbL5cXksMZ7QR21T+Pno0av78M9t1I+rKFfAnro2KQoZp05DriSdMCxQO9La2+rPOiG0FOBuOu+c8BUjKIImEQYD0DZuq2Vy+fNl8+kmbNi2KFy+uACkAS2zKIEm442tLhgwZcPHixYAGSLRrzx7UfO45/JwpE/INGoQyADKwPqlvX6T97js0LlUKszdsMJktDt6eNGkS+vfvbwrDd+7caWqbWGjOAdsWBUgiYRAgrVu37qrruOT00EMPoU2bNr44LkmKltgkAnHHWZkyZcyyFQMPBkoMmGjNmjV+D5BuueWWq5bDfoqJwRsAXgPQ17px0SJMy5ABXdk/yeHA8uXL8fLLL7tqpXj8DJDs9UfWY/J2tU8RCQ4++0vkGzSbRyp7lAq0xCYRvszGTRJbt2415xmA8OSrAu1kZZCcARIzVxfYnZuz2D75hAdgrs9y8SJmAJgI4KkePRIspR05csR8tV9HzChxnpuIBAefflRhZ+3EBjmKD2mJTSKUe6E2u1CzMNpqsti9e3e//WxPAZI12iRHjhzI36kTt/gCtuLxRwHM2LULVWyPc/jwYY8BEmmZTSTEl9jGjRuX4DILEPnJZ+rUqWjevLmvjk0SowySRCh7gMQC6IkTJ2L//v3mcsOGDc0YJH/Jnj27ObH7NgMk/tyDBw+a26pUqRI/kzJHDrbHxvn69RHbsyeyArgVwAoAzwF4M4kMkhUg1a1b12/PQUT8HCCNHj06wWWumbP3BztsDxkyxJuHlJRQDZJEKHa4towdO9bUIVHBggVNN2tuFPEnZpFYg7l7926sWrXKc91TVBQyPfkkGr/2Gl7bswe3A+CQlJEA7gTw7b//mru51yDZM0hcQmQAZe+plBQGa/fccw9y5cqF2bNnu2qzRMR7Xr2aRHK/jqBoFKkMkkQoFjdnzpwZ586dcwVHDIo4ssNftUeeAiT+/c+aNSvJwvBs1aqh1p49eAnAYGc9QxMA1Tmq5NtvPWaQ+NrK5UJ25GZ905QpU9C5c+drHhez91bPpblz56JVq1bX/VxFIl2Ka5C4pb9x48am1X8kYpNIjhFYvXp10NQgseeK+7gXkXDEbDV3srmPNqpVq1aq/Hx7HZI9QOISm7t27dqZ7trjCxTAyZkzsdd5/Y2XLwNt2qD21KnI7CGD9Msvv5hdeSxd+PTTT5N1XNYyI11r/puI+ClA4hvxxo0bU/pt4sclNmWPJJJUqFDBdZ6jPXr16pVqP9seIFlLZOyozx5N7nhsW7ZsMcXkOdu2RbP8+c3ONsvta9eCzQlucwuQfvjhhxRn663Cb7LqokQkALvYOnXqhA8//PA6f7T4aolN9UcSSdhwkYFKixYtTJG2KY4OQIBk4VzKxI6hVKlSyJkzpzmfLl8+tAPwWJo0cGSOzx0xF/YbgOezZAEfgZtdZs6c6fr+f/75J1nL+QqQRIKkBolr5JMnT8bChQvN2nuWLGZ0Y4KUt6TeElshZZAkglSsWDFg2+E9BUjJbUyZN29e83VibCzeXLIEh5o0QclTp0wB97CzZ8EBJA85gyJ7SQOXz4oUKZLkY7PdgUUBkkgAM0hMGXPNnUs727ZtM0WL1ulawxnFdxmkOADnlUESSdUicW8DJPuOtEPZs+OhW2/F67bb7+VwWwDum/ytYHDRokWmjcp3332XZAbJHiyJSCpnkBYvXnwdP1J8FSCdYw8q1SCJpJpMmTKZgbP2LI2nAu1rBUjshXT45Emzu21V1qz4NDoamU6dQmG+vgJ4EcBw54cgBkj169c3tVasaWINaMuWLV2PxSW4o0ePui4rgyQSwAwSe4Bwh0Vit0nqLLGpi7ZIYJfZ+OGkZMmSKQ6QmPGxirzX582LDR9/jEXO29IA4OS2+QDyOwu1OfCWwRFxyc3eQ+nYsWMJXo8ZICX2+iwifg6Q+AJhdYO14x+qpzV68U8GSU0iRVKffWgtC7STO1zWqkFyD5BYxF349ttNE0lOsrRKshuzWzjnuv3yi2ktYsdhvfbHcp/pppFPIgEKkPjpxNOuDX7KYU8eSZ0ASU0iRVKf/UNgcpfX3DNIDHDi4uJcARI7gadNnx6vcGQKawudO98YUj29eDEyvPhignqIpAIk0jKbSCrXIA0YMMB8ZXD0/PPPm4629nXw3377DbfdZu/qEX4C3kmbP/c8S7O1xCYSCPbXuHr16nkVIG3dutV1ngESs1C33nqr2QDzO19X16/HgpIl0dTZLbzSjz/iZ/ZWYhlDMgOk0qVLe/X8RMSLAIm71KwM0qZNm5A+PTeoxuP5SpUqYdCgQQj3Tto8nTp1yjSIS3XnWJodTxkkkdTXunVrvPLKK2bESUpGetiX2NwDJHr99dfN4/bo0QOZihTBS1Wq4McVK8xON/bJr+Xc5fawMkgiwRcgWbvXHn74YYwbN05vzIGgQbUiAZUmTRo8++yzKf4+ewZp586drvM5cuQwX++++25zstxcvDhGr1iBXwB8ydon3pfDbgF8NX8+wOxShgweA6TkbPVnrdKIESPMh90333zTYwsDkUjmVQ0Sd21wOKQ7No/kpyDxIw2qFQlJzDhbMxPZbNc9g5RYrROnPrLSyf6K244ZIs6f27nTYzB0rQzS3r17zfLgiy++aDp3Dxs2zMtnJRK+vAqQJkyY4HF9u1y5cnj//fd9cVySzC7apAySSPBj7Wbu3Lmvuj6xAMm+W4570u4H0APABevKtWvZpRLFNnCvW/IDpLNnz6J27dqmZtSybNmylD0ZkQjgVYDEP74CBQp4TCFzlpCk7hKbMkgiocFeh5TcDJLdxxkzogZrmKwrYmLw7OrVpk6J/ZOSEyDNnz//qn51rGmyN5sUES8DJM4FWr58+VXX8zpuV5XUXWJTBkkkNNjrkFKSQbLceeed2AjgdgB77rjDdf3/APycJg0Kp0lzzRokduK23HTTTa7zK1euTMEzEQl/XgVI3bt3R79+/fDRRx/h33//NSfWH3HKNm+T1F1iUwZJJHQDJKtI2x0/bNp3CpO1a+40gKn33guMGYPLzttqx8bid4cD9a+RQWJRtuWxxx5znVeAJOKDAOmpp55Ct27d8OSTT5pPOTz17t0bffr0wZAhQ7x5SEku7WITiYgMEnfL2TM83GXGzt2WHTt34lz37mAnpj3O6/LGxZmRJQ8dPIg4WyG4pwwS58p16tTJdf2KFSu8fl4i4Sja22JD7lbjuBF+6tiwYQOOHz+OF154wfdHKAlpF5tIWNUgJZZBcl9mK1++PEqUKJGgbsi8Bjt3uW1wPjYX2YY7HLjMDJNtZhudO3fO1UOJm2qKFi2KwoU5IhdYtWpV4BrgioRLgGTJmjUrbr/9dvOHmyFDBt8dlSQ7QOKnQDasE5HQyyBxGkFS45nshdoVKlQwH4by5ctnLjPQsXogsbx6Ytu2mHXbbYgfYAJkmDePs1CA3393Pcaff/7pGmTLx6M7nLVMHBXF20XkOgMkbgtlerZWrVrYt2+fuW7q1Kn45Re2NZPUqkHS8ppI6AZIiS2vWUqVKnXViBMri8Qdw3///bfr9rwFCmDV3XejuTNgMv75J75f0oQJHIGQoP6oYsWKCQIkT3VI33zzjVkZsAbrikQSrwIkNhZr1qyZyV6sXbsWF53zgjhBevjw4b4+RkmiBknLayKhGyAltbxGXbt2RfPmzdGhQwe0adPGXGdfZrPXDXH5jtml+c4ltyPFi8ffcOkS8PjjwEMP4S9bNsnKINWsWdNjgMRxKG3btsXLL79sOm2LRBqvAiTOC2JDyIkTJ7o6wxKbjzFgCmccVFu2bFmztBgQbktsyiCJhG4N0rUySAyg5syZgy+++ML1WmsPkOztVvjY+fPnN+dZtD310UeBPn3+e7BPPsFjn3yCm9wySKxFstj7I3344Yeu8/rgK5HIqwCJnyw8TbFmK/2TJ08inHFQ7ebNm7F6NQcABH6JTRkkkfBdYrvWstuaNWs8Bki0n40fx44FvvwSyJLFXFfyzBnwO+7Lnt11LPyQZbUTsM91sz+2akwlEnkVIPGP0D5N2sL6I0/NzcR/S2zKIImEjuzZsyfYVOFNgMQlNzbrdeceIHG4+OXLl4H27XFs9mycd+5W47CTL/lBlstmDofZlWwFS9wVR6dPn05QT8oNOSKRxutGkX379jWzfPjHtX//fnz22WcYNGgQevTgtCBJrSU2ZZBEQoc9GElODZInDFY++OADjwESl9+sOZksd+Aw2oEDByJPw4YosHcvZjvva/ptP/UU0LGjeU2xlv44biQuLg4LFy6MD66cjh07hgsXLrgG7X777bcJOnKLhCOv9ocPHjzY/BE1btzY9NXgchtTsAyQ2DBS/Ei72ERCmn1mpTcZJCuLxALuKVOmmMtcImOJAwMwXsfdxexpZK8d4sDbFgBeBODqWDdtGrB5Mypmy4Z1zuCHZRKse3LH3crFixc3wVmvXr1Mi4Lt27drvJSELa8bRT777LOmOeQff/xhdj4wNcvdDpI6GSR+trukDJJIyLFnkLwNkGj06NGuJTXWJfF1mapXr24yR3Zc1uvYsSPatG2LuKFDETdzJl884m/ctAnvrFqFZs77sg7JU4C0Z88e17Bb4ofjH3/80evjFwl219VhkJ9aypQpY85bf5ySOgGSBtWKRHaAxOW5RYsWmd1mXbp0SXAbRz799NNPpg6JW/9nzJiBOnXqJHyAsmWB1q256wZZLl0CQ6LnAPy+Zo0pm3C3d+9e89XeTJKPz7FTIuHI60aR/KNkB212geWJ5ydNmuTbo5NrBkjKIImE7lb/6wmQiC1H3nrrLVSqVOmqOW7M9CxYsMDsOr4qOCLWKq1axQm4rjcDLshVeOklWCXZhQoVShAgnT17NkFzSgZhVmdukXDjVYDEzqos0m7RogW++uorc+L5/v37ax5bKtUgKYMkEpratWtnajY5RoS1Qv7CZbUmTZqY2qRE8fXj66+xplUr14iSijt24FcOx3VrIskA6a+//koQELGWatu2bX57DiIht8T23nvvmSaRXNO2tGzZ0jQeY5H2sGHDfHmM4iGDZJVqK4MkElqYzTl48CCyZMmSoNFuwERHY/8jj2DorFn4jK0ImEUCwE5vP+XNixm2AIk1p+6YRbL3ZhKJ6AwSt39Wq1btquurVq1qdkGIn3BkgHPrrTJIIqHdDykogiNbXRRrkKqzEbDzOvZLuu+DD/CEs740sQCJdUgi4cirAKlz584mi+RuwoQJePDBB31xXJKMHkikAElEfFUXtR1ADQBznddHx8biPYcDb3M5bffuBAGS1V2bAZLqkCQcpb2eIm0WAVqToNk0knN8uJtiwIABrvuNGjXKN0cqV3XRJi2xiYgvd9axX9K9AN4A0N95XS8AZY4cQQ/n+BHuoGNJxc8//2yaS7J4W922Jdx4FSDxU0SVKpwXDezcudN8zZ07tznZP2Fo67+PKYMkIn7AD1ps23KJy/gAYgG8VagQ+g8bhiuPPoq0DgcaA5hz7Bha8vW+fPkEu/EYJClAknDjVYCkNecAceuiTcogicj14odZBjxWryMqzNltjzyCd+fNQ/vp05EPQAkAKwFMveEG/JmbVUr/BUjFinHfG8e7xc93E4nIGqTz58+bLqqWf//9F2PGjHF1WA1n77zzjuk9cvvtt6f+D9cSm4ikwjIbWQNxr9SoAb7acRQJMWf95I8/opVte7815HbdunW46aabTNsXjqMSibgAqVWrVvjkk0/Mec7tYWt7Nivj9Z6Kt8NJz549sXnzZqxezU2wgV1i4ywkNoQTEfF1gGQySM6vHDLCVpPWlv8ohwPNFi3CROcyBDNI1kYdjiT54YcfsMZZryQSUQESp0TXrVvXnGcLe84DYhaJQdO4ceN8fYySSICk+iMR8RV7TZE9g8RmkWxJwDWDj+66C2cHDXLd51EAnMZ2avduc5kbdSyexpWIhH0NEpfXrKUdLqu1bdsW0dHRZkcbAyXxfw0Sz2l5TUT8nUFioMSls0OHDqFBgwbmtR5VqiC2a1ekuXwZTdh9e/Ro4IEHEtQw8f4iEZdBKlGiBL799luTSp03bx7uvPNO1xRoZTX8SBkkEUnlGiQqV64cGjVqFB8cUceO2DlhAuIrj4C8x44BNWog/65dru9ht3CRiJzFNmjQILNroUaNGq55PcwmVa5c2dfHKIkESMogiYi/l9gSk6VpU7AL3l/WFUeOYNbp07jfeVEZJInIJbb77rvPzBPioEL7FOnGjRujTZs2vjw+SWKbvzJIIuKPDBIzRawtTUquXLnwNwCO212YLRuqxsQgI4BpAIpzZIkzg8TdbK7Mk0gI8fq3ln88zBbZf/G5m6106dK+Oja5xjZ/ZZBExB8BUsGCBZE2bdKfnzNmzGiaQ54E8FC+fDjQvLnrtuG8bulS1LvjDhQqVMjjDDeRYKewPpSoBklEUmGJzSrQTm5QdeDYMSzq2BGDbbe1OHYML/32Gy4ePIgpU6b4/HhF/E0BUggvsSmDJCK+wqAoX758rq39ycHxUnT8+HH8u3s3XgfQjs2Enbc3BLACwKW/XJVKIuE/rFYCv8SmDJKI+EqGDBmwaNEirFy5Evffb5VaJy+DxPEiGzduNOfZTJLdkL4DzHiSUgCGcsoCm+sGYgKBSGplkC5fvmyKsbdv3+7tzxRvaYlNRPyI2/m7deuW7Oy0lUGi9evXu86vAlADwJ/OyzkvXwYaNADmzPH4OGfPnkWHDh3QpUsX18BckZALkNhR1fqkIKlMS2wiEqSF3dtss9noX+d4kp+tKzi/s2VL4MMPr3qcl19+GdOmTcPUqVMxefJkfx+2iP9qkDp16oQPPfySS+pkkLi+zzGQyiCJSCDZM0iecIdbMwDTrStiY4FHHwVeeonrcq772d9P2HxYJGRrkK5cuWKi/IULF6Jq1arIkiVLgttHjRrlq+MTDwGStdCmDJKIBHOARBcBdADQqFMn5P700/grX3wR4FgSDjdPmxanTp1y3T9btmz+PGQR/wZI7GlRpUoVj2nVqKgobx5SvAiQlEESkWAaT0KZM2c28zrtmCta27kz7uT7xoAB8VdOmsSJtsC0aQnqjmKZZRIJ1QBp8eLFvj8SSXYNkjJIIhKsGaTbbrsNv/7661XX79u3D+jfHyhUCOjcGWBQNGcO4urXB8Msa66bRpRIyPdBWrZsmalFqlWrVvwvPmAK7H755RdfHp9YuF7vzCBZpdrKIIlIsGWQ7r33Xo/33c9sEbGFALf9O5fSoteuNb2SSjjvpyG3EtIB0syZM9GsWTNkypQJa9euxcWLXGUGYmJiMHw4m8yLz1244CpqVAZJRIIxg9SiRQv069fP432tD9JG/frA8uXsTmkucnYbc07VlUGSUA+QXnnlFbz//vuYOHGi2fZvqV27tgmYwtk777yDsmXL4vbUbnjmtsWfOAdJRCRQsmfPjltuucXVQ+mzzz4zH5yTzCA5eyYN+ugjbJ0yBYed3buZi2Lxxu2HD5uNQM8884yZ7/n777+n0rMR8UGAtHXrVtSrV++q67n74ORJbuwMXz179sTmzZuxml1hA9hFmzsH06RJk7rHICJiw2Hls2fPNh8cWV5hZbWtLFL//v1dA82tDNLOnTtRt25dvPXWW+j6zDMY2qgRFjkfLzOAb9lgcvBgjBgxwrzODrCKukVCIUDKnz8/duzYcdX1/AOxPk2Ij6mLtogEodKlS+PJJ5802STL6NGjcezYMdPyhe8XVoDE3WodO3bEGWdG/LfffsPqbdvQHMCXtp1Dld56CwOdl5cuXWpGmYiERIDUvXt39O3b1/xyc1s/U6dMrQ4aNAg9evTw/VGKumiLSEjJmTOn+VqwYEHz9cCBA2jfvv1V2XcuoV0G8ACAt23XvwmY4bee2smIBO02/8GDByMuLs7MZGO/Cy63cdAhA6TevXv7/ihFg2pFJCQVKlQIa9asMee//ZYLaJ4xR9SHRdqsc3Ve9z8AeQGsWLoUpUqVcg3GZf3r6dOnzVKeSg0kqDJIzBo9++yzOH78uGkayenPR44cMfN0xE+0xCYiIYibd+z4YZqrEO64+YVeBfAYG0Y6r38IQBW+tzibT7KU4/HHHzcfyFn7JBJ0fZAoffr05peaOw20oyp1AyQtsYlIKBg4cCCmT59u+ubVr18fc+bM8dgKoHlzViLFmwignXNMCVXcswdo2hQ4fhzz2UPJiaUeqk+SoAqQunTpYmaxcTeCBKYGSRkkEQkF3MXWrl0700h4yZIlaNSoEUqUKIG0aRNWeNx9990JLp9q3Bj9SpWCa0obu3PXq4coez8lAKtWrfL3U5AIFe1t5ui1115DyZIlUaRIEfPJYNKkSdi+fbvvj1A81iApgyQioYrvIQySLDly5EC1atUS3KdXr16IbtwY9dld27ryzz/R68svEV+NFO+jjz5KnYOWiONVgMRgiLsK9uzZg5EjR5rlNfa04HbPws7OqOJjqkESkTDCvnn2AIkf+nLlypVgZAnnuq1nHRMnNThvy3v+PDjQymrVy93UIkFXg8Rfav5C8yt7YDBl6mk2j/iAtvmLSJhlkSx8H+HmH05o4FIbZ33y/YQBEv0NYGiTJnA4L3PAyU8AmgL4999/A/YcJLx5FSCxBTyH1PKXmlv+L1y4YL5yyOC6det8f5Sibf4iElZeeOEF13mWbNB9991nOnPXqVPHXC5fvryrE/fSbdtw6MsvscT5PdwW9AOAO0+cwKlTrkol08SYEw/mzZuXqs9Hwo9XfZD4y8xM0dChQ9G2bVvceuutvj8ySUi72EQkjLCP3tdff20CIBZue8K5bizd4HgnfvguW7MmuNn/cwBtmYVynj/42mu40Tko/amnnjL9lt59913TfsZ9oK6IXzNI/EVlHyTuHmCPCzYCe+CBBzBhwgR1PPUX7WITkTDCJbU2bdqgVatWSd7PWmajEydOmK3/bAHwsXNQOt/ECo4YAQwdyi6SCZpRcuecSKoGSJUqVUKfPn1M9M8InX0tuJ7MtGaZMmW8PhhJgnaxiUgEsgdIljh25W7e3DSVdBk2DHFPPpngTe3DDz9UnyRJ3QCJv3Br1641gwhbtmyJhg0b4tNPP0WFChVM4CR+oF1sIhKB7rnnHo/jRBo1boznANhbTka//75ZcrPKv//8808z7UEk1QIkDiGsUaMGPv/8c9MLacqUKTh69KgJmjjFWfwXIPGT03llkEQkQnBaw8aNG81QW77XcLWCRdxWY8mxAMbdfjsczsaT7Z3F21mc3//PP/8E8Ogl4oq0mS2qW7eushgBqEGy8kj6txeRSGHNaatSpYqpW2LvvcuXL5sCbw5Of+/0aazNmRPvHj6MzM7t/2wDwBCKZSB09uxZZMlihU0ifsogMeXJX0o2h3z00UfNicttMTEx3jycpCCDpABJRCIZs+cs8GYmqWDBgua6LVu2YMrhwyYwOuG8X3UAyzjPbft2PPTQQ+Y103247cKFC9UOQBIV5fCigm3NmjVo1qyZ2YLJQbW0evVqnD9/3gwSZJQf7th3g51gGRSmSrDCBpxHj4LT79igPzY21tUfREQkEnEl45df2Ff7P+UBMOSJD52AE1mzotaZM9jivGy95bED9x133GHOf/PNN2jdunWqHrsE//u3V++w/fv3N8XZXNvlTjaedu3aZVrDe5rSLL5dYuMnKAVHIhLp2GbGHUuyawGwJoPmOHPGZJKs0SRc/aDvvvvO9T1cthNx59W7LDNITz/9dIJpzDz/v//9z9wmPhYbC1y4YM5qi7+IyH/lHp5w+AhDp7XOy9ZokiYA9u7da65zXzzZunWr349XIiBAYkpq9+7dV13P4bV68/aDc+wdG09b/EVE4tWsWTPR649GRaGBMzCyRpPM5tDbCRPM5UOHDiX4npkzZ/r9eCUCAqT27dujW7dumDZtmgmKePryyy9NsXbHjh19f5SRToNqRUSuwpWLW265xXW5XLlyaNeuHcaOHWva0Zx27mSzQh/2RyrPkSTvvntVgOTpQ79ENq+2+b/55ptmF0GXLl1w5coVc126dOnQo0cP19BB8SENqhUR8YjdstmsOGPGjPjhhx9QrFgxcz3nhR47dsyMJrkfwPsAusfvTAJ69kTLAgVMRsly4MCBgD0HCaMAidsrGaGPGDECO3dyXxVQvHhxZM7MDhTicxpUKyLiUYMGDbB9+3bT46hAgQKu6xkgcfs/sSz7MfZEAvCM8/bHDhzAJQCc/cBqJAVI4pMAycKAqHz58q7Bg+InGjMiIpKoEiXY/CQhBkjunnUGSda8h17OAu4uCpDEg+jrSWsyOGJakyeenzRpEsIdG42xq+vtt1ubRlO3Bkm72ERErs1TgERjnAFRfHEI0AHA9wBOHzjgagEg4nWA9MILL6Bv375o0aIFvvrqK3PiefZH4m3hrGfPnti8ebNpjJlqlEESEUmR3LmZG/JsKoDWzrmW1AzA3NhYHNu2DcuXL8dzzz2Hffv2pdqxShgtsb333nuYOHFigh1rbBxZsWJF9O7dG8OGDfPlMYpbgFRIGSQRkRRnkF599VU8+ywX2uK3/Dd1Zo9yAGBP7fN3340nT57ExhMnsHLlSjOKRCKXVxkkDgmsVq3aVddXrVrVtatN/LfNXxkkEZGk5c+fP8FlliYUKlQowXXLAdyVKRP2Oy9n2rUL3584gVIAFi1aZIIkvadFLq8CpM6dO5sskrsJEybgwQcf9MVxiZ22+YuIpIg1Z82SK1cuj1ml9FWrmq7b1miSogA43e12Z8NJToiQyHTdRdpsDslThQoVzLIbZ4QNGDDAdRIf0DZ/EZEUuemmm0z7GQs3E+XNm/eq+9WoUQP/AKjDcSNZslw1mmTcuHFmOPjFixfx119/XTWiRMKXVwHSH3/8gSpVqphonH2QeGJBHK/jbevWrTOn9evX+/6II5GW2EREUsy+25g9kTxlkOrUYWgEHGZPpbg4LHEbTdI2NhabNm1C8+bNzQ7mIUOGpNrxSwgWaS9evNj3RyLJXmJTBklE5Npq1aplxmARR4+4B0hsLsnxJJaD58/jLgBfAGjjHE3C7576+ONYvGqVuc/rr7/ucWLEpUuXzI5ujj5JbEZcYpiV4qD3IkWKXFU7JSG4xCapSNv8RURS7KGHHjJBB40ePdo0N86QIYPrdi7BsTbJjqNJ2gH4q1Yt15tk11Wr8PI1ftYHH3yATp06maBs/36r7Dt5GMRVr14dpUuXxqlTp1L0veI/CpBCgWqQRERSjK+VGzduNL3rODuUWEtkD5CyZct21ffFMns0ZQpGZ8zouu45AJOcyy4XLly46nv69OHQknjTp09P0XE+8MAD5mtMTAw+//zzFH2v+I8CpFDg1klbGSQRkeTJnj07ypQp4/G2m2++GWnSpDH3cVegYEH83KyZGUdi9dfuBuBbAJ3atDH1Taucy27u0qb1fooX2+hIcFCAFGIZpHPOdXMREbk+OXLkcNUnuWeeuBxXv359vAPgfufSG90D4Km5c/HPmjUYPny4x8dNly6d349dgjBAYnTbuHFjMz1ZUjdA4ueKjDfeqMHAIiI+YG37dw+QrELpBg0amK8zneNITjpvr+FsMnk0kZFTzEpZzpw5Y2qh+vXrp1lv4R4gMTLmmq6k/hKb6o9ERK7Pp59+ar6yOLtDhw5JBkgcn2Wd/xlAXQDWhLZbAczYvx+XVq686mecP29NeQOGDh2KKVOmYOzYsZg2bdo1j099lkJ8iY2V+mwUKambQVL9kYjI9WFB9Nq1a03htvV66h4gFSxY0JUJmjFjhinAHjNmDP5gd20Am533Y+iUpnFjYMGCBN9/1lYW8e6777rOJ2e2mwKk4OFVJRln00yePNn8z+b8NfeamFGjRvnq+IR/MGfPgotqyiCJiFwflihUrlzZYy2Sxd6Bu3bt2uZkrZzscXbd/s75Nc25c4hr3hzch/a5bVnNYt/xlpz6UfsuOwnBAMnqpE3btm1LcJvqY/y7xKYMkoiIb7lnkNjs0V2BAgVc508AaOoMiNhQMjo2Fp/xPgDecssg2bHw+1rOneNWHAkG6qQd7C5dQpRzmrSW2EREAhMgsWaJTSf37Nljum//+eefuA/A2wCedN7nTQCFONbElkFKrDYpMQqQgoe2+Qc7NYkUEQl4gMRB7HPmzDElJD/99BMefvhh0x+pJ4Bnbffrzw7e8+dzreyqXWsnT1r74BKnACl4eN/NioVqmzdj9+7dZgaNXcuWLa/3uMSiMSMiIqkaIBUuXNjj/cqXL29OxDrcokWL4qWXXgK7IR3kuBHnm2rN3btxsEoVPOzWpVsBUgQESH///TfatGljJhyz5siqurfqj2Jj2ahd/NFFWxkkERHfci/StvcxSgr7G3F4LQuxJzuDpK9Ya8Qdbps3gyNtNwA4YKvfZWDFJELu3Lk97lpTgBTiS2x9+/Y1LdoPHz5sis64Frt06VJUq1YNS5Ys8f1RRjJlkERE/Mrb0SDFihXDX3/95bo8B0BDNpB0Xq4E4FcApWzJhW7dupnAyuI+100BUogHSCtWrMCwYcNMBMx1WZ7q1KmDESNGJBjYJz6gGiQREb/irmyrIeRbb3EfWsqCJDtOZ6sFYJd1u7Pr9h22+8yePRtHjx71GBApQArxAIlLaNYbNYOk/fv3m/M33XQTtm7d6tsjjHS2JTZlkEREfC9DhgxYt26d2aHNFZKUGjhwYILL250NJdc6L+cCsAjA3bb7sOD7vffeS9A+ILk73SR1eJVXZJHahg0bzDJbjRo1MHLkSKRPnx4TJkzwWP0vvskgqQZJRMQ/mEGyskgpxRWVU6dOYeLEia7rDnGWm3OOW1NnXdIsAF2d/ZO4E44jSNwpgxTiGaTnnnvOtX2Rvxi7du1C3bp1TUQ8btw4Xx9jZFMNkohIUGMtLhME7k4DKLdrF352Bl7MSLChZC8WdB9kSffVFCCFeAapWTPONY5XokQJbNmyBcePHzc7AdRJ279LbMogiYiEjgI33YRB9eph8/Tp6OG8js0lX5w3D/M83F8BUhg2imQfCQVH/l9iUwZJRCQ0tG/f3rwvFipa1HTbftl224sAuN7i/q6ZnF5JEmQZpAEDBiT7QTWs1oe0i01EJOQ8//zz5kT16tXDm2++iRcAHAMwxnmf3s4CbtYlXbENuv3222/RunXrgB27pDBAYoW/BIB2sYmIhIQmTZpg4cKF5nz37t2RLl06c541upaxAI6zE7fzDfgBANkBM9fN2r/G2t67777bbH6SEAiQNKA28Bmkc1FRyZoGLSIiqY+72IYOHYoGDRqYwbaW7Nmzm2DHGss1lUtpAKYDyOjc/j8fwL0AYpwJCdb3sgmzVg1CrEib0W1iuN5qpRXFtwESsmRRnZeISJBi00hPW/dp2bJl5r1z3rx5uHLlCr7nhicA3wHgxLY67LqdNi0aX7liRpbs2bPH9En63//+l+rPQ+JFOdwHwSRD5cqVE1y+fPmy2erPdu3FixfH2rVWe6zwxZ4X2bJlQ0xMjH+Xvdq0Ab791pytVrAg1uzb57+fJSIiftWrVy+88847rsu3AWY3W17n5Z3OvknsxF29enUcOnQI5cqVw/fff2+mVkjqvX97lUHyVI/EH8j5MhxiK/6pQYpSqlVEJKQxe2S3HsCOjz9G3hdeAHbvRnHnaJI7ObZkFQeXAP/++y+WL1+eoJZJ/M9n4SijsJdeeknLaz7msAVI6bOzlE9ERELVxYsXr7ouY4UKwPLl2OPMZnD4yFLnuBILd8L17t3bDIWfP58VS+JvPs3XMV3Fk/hO3Gn2Yo3f3ZBFO9hEREKap3FcZpmncGGMadMGvzmvywGA++H+a8sMjB8/Hr///ruZZkFHjhwxNU8HDhxIpaOPLF4tsbmPE2EZE/8HTZ06Fc2bN/fVsQn/bZ0Bkrb4i4iEvn79+uGzzz5LMNg9X7585mu6/PnRGMA3tvltLObuDGCa7TFWr15tvj744INYsGCBySpxOU6beIIgQBo9enSCyywcy5MnD7p27YohQ4b46tjEtotNg2pFREIfX8c3b96MHTt2mEHvTCpYr+0McM46t/t/CqAdgybncFtmlN63PQ6LjFn7S2vWrMEff/yBClyqk8AGSNyxJqkj6nx86zBlkEREwgOTCrfeeismTZrksYCb3ZI6AMhapAia79ljamHeA5AbwCvO+1rBkWX69OkKkHxMewaDmcOBaFuApAySiEj44kqMpdgtt6D5v/9ihO12znLj+o2nhbSvvvrKlLtIgDNIic1lY3owY8aMpgNoq1atzABbuQ7nz7NRlTmrQbUiIuHtiSeeMCUsx48fx+eff843VWR46y0MGjgQbzrv04/D4QF0s81vI9Y0MTPFrt0sdeGOt0yZMgXomURwo8iGDRuaZpCxsbEoVaqUuW7btm1IkyYNSpcubf5HMVj65ZdfULZsWQQrDg/86KOPzLEOHjwYnTp1Cq5GkUeOAHnj24exUG/vu++iR48e/vlZIiIScBcuXMD58+eRIwerjuJVrVoVFdeuBRfk0jivYwfu9rx/Io/DHW89e/Z0ZZVUwJ3y92+vltiYHeJQvv3795sthzzt3bsXTZs2RceOHbFv3z7Ts6F///4IVps2bTIROo+dOwL4y3TyJKfjBOeYEdUgiYiEP67C2IMjy8fOgbZWF6WWzg7cfFfwdP8vvvgCO3fuNDPhKlWqhLP2sVWSLF4FSG+88QZefvnlBG/YjMZefPFFU5XPgaovvPCCCT6C1V9//YWaNWuaX0amIfkLNHfuXAQVW5NI7WITEYlM6dJxLxvAoVN3AYhv/gLUA/AzgCqFCl31PXnz5kWHDh1MwoIJAQZMkgoBEtNShw8fvup6Nq2yKuu5DmpNLvbG0qVL0aJFCxQsWNCkBr91ziOz4zwbDgdkkFOjRg1XW/bkKF++PJYsWWKyRidOnDDn+YsUVJRBEhGJeCwHsTwzfz4+6twZR2yz3D7asQM3eXif5vZ/i6f3bPFDkTaX2B555BG89dZbuP322811XKYaNGgQWrdubS4zWOE2Rm8xHcisDn9O27Ztr7p92rRpplj8/fffN8HRmDFj0KxZM1P/xMiZbrvttqvm3hDbtLM2qk+fPmjUqJHJft1xxx2mhiqoKEASEYl4derUwc8//2wySVz5YM+julOnggNHigIocuGCySQ1dA65JS6vuQ+VlxRyeOH06dOORx991JE+fXpHdHS0OfF89+7dHWfOnDH3WbdunTn5Ag/zm2++SXBd9erVHT179nRdjo2NdRQsWNAxYsQIr35Gt27dHD/88EOit1+4cMERExPjOu3Zs8ccF8/7zaxZfPLm9Azg2LZtm/9+loiIhIRly5aZ95/CgGNv1qyu94ndgONmwNzmfurRo0egDzto8H07Oe/fXi2xZc2aFRMnTsSxY8ewbt06c+L5CRMmIEuWLK7sDU/+wKU71jexUNzC7Y28vGLFimQ/jpVyZNaJGS9moBIzYsQIk2myTix8S80Mkrb5i4iIlVFine8d992HtL/8gn3O9wa+K/3obAPgTktsqbTEZg+UKlasiNR29OhR02LAml9j4eUtW7akaKmQ67QM6rjdP23axP852FfC3v+JtVZ+D5LclthUpC0iIvTSSy+5zg9o0QKPfvYZ2FSnlLMFQBO3FgAzZ87E7t27UbQoF+XErwHSokWLzIlRaVxcXILbJk+ejFCQkmxThgwZzClQu9jOR0Wp6ZeIiFwlNlcucEz8SgAFANQGMNXZJ8n+7vzhhx8mCKwkaV4tsfEf+M477zQBErM53AVmP/lb7ty5TUH1oUOHElzPy/nz50fYsGWQ4jJnVqMvERG5Clvr7AZwj7McA86eSf/tfYu3YcOGABxdhGWQuHPs448/RufOnREI6dOnN51FGaBZu+aYxeLlXr16IWzYG3s5a7tERETcAyRa5wyMZkdFIY3DAbZq/hfAWNt7p/g5g8Qi6Vq1asGfzpw5g/Xr15sT7dq1y5znGiqxHoiF4lOmTDFNHzmCg60BHn74YYQNBUgiIpLMAAnO7trLbWOzxgB40HnemhbBzeHsWyh+CJAeffTR+EF6fsQGV5UrVzYnKyDieVbuU/v27U3zLF7mbjkGT+yE7V64HcriTlv9UoFo7WATEREP3DNDuZ5+GsPcxpS0tAVI99xzj3mvZGKB5ydN4pS3/zARccVDD8FIk9bbYXrc0r9w4UKzi81qg24ZNWrUdR9YgwYNXEP2EsPltLBaUnNz5eRJWL/2aRQgiYhIMpQrVw6VoqORKy4OPZ1v9NMBPLJnj2mU/OOPP7rKZWjOnDl44IEHTCbqgw8+wBNPPGEaMHMjUyTXvnoVIG3cuNHV44gdPe0i+R/T12KdY1soXfbsAT0WEREJHWt+/x1PP/UUmsTEoNTq1eAe7A8OHkTjRPr9sXUNAyQGR/Tbb79hx44dKFmyJCKVVwHS4sWLfX8kcpU4BUgiIuIFJjHmLVjAGSNYkicPGsTEIKuzkSSH3P7poe7X3cWLFxHJrqtR5ObNm81apX0oLTNIHDIr18/h/IVlH4uMOXIE+nBERCTUpEuHEZUqIW7pUjRydtleCJjzf9nupgDJRwHS33//jTZt2mDTpk0mILJqhazlNXa5Ft/tYjODarNlC/TRiIhIEOJgd0vDhhxZm1CWXLnQyhkY1QDAboE/uQVJ27dvNyO7cI2gKZJ4tYutb9++uPnmm00Xba5Z/vnnn1i6dCmqVauGJUuW+P4oI1TUuXPmq8aMiIhIYurVq2fGYd17772mR6G7PHnymAaSd7E2yXmdFSSVcV6+//77EwRadNq2kzoSeZVBYmX7Tz/9ZDpaM+LkicPzONC1T58+ZnitXL80F+In6WhQrYiIJGX48OGJ3la+fHnzlZv8mwJYAKBaIpkku7179+KRRx5B3rx5za42ttbp168f2rVrh0jgVYDEJTQro8Egaf/+/ShVqhRuuukmbN261dfHGLHSOtd/lUESERFvWbvOEwuSFjuDpM0eVosu2WqM6ddff71mC56IXmJjNGrNdGFUOXLkSCxfvhzDhg3DLbfc4utjjEyxsUh7+fJ/NUjKIImIiBfYr9DOCpKs5Ta2V17qDJjsLrkFR5Zjx44hEngVID333HNm9hkxKOIYkLp165pmU+PGjfP1MUYm25gRLrEpgyQiIt7Ili2bSWS4B0lfP/EEVjkv53IutzVIxuP99ZenBbnw41WA1KxZM7Rt29acL1GiBLZs2YKjR4+aou1GjZioE18GSMogiYjI9XjqqadMax677DffjCbOJTa6wdkniWNJksL3/EjgVYDkSc6cOdVF248BkjJIIiJyPcqUsfasxStUqBC4T+1uAN85r8sIYCaA/8bdXs0aGh/ufBYgiY/Z+k9oF5uIiPgCB7tXrVrV7DgvXry4uY77pf8PwKe23VtTAfRP5DFefvnlsJ6DaolyREo5uo9xbg3XdWNiYvwTvCxfDtSpY86+BeDJc+eQKVMm3/8cERGJSHz7Z2PJn3/+2VzmGhCriO2hzzhnoBRfdZwQd7AXKFAA4fr+rQxSCCyxnY+KQsaMTHyKiIj4BstiOFv14MGD6NChA3Llzo3eAIba7tPHueTm6eP57NmzEc4UIIXAEtuVjBlV3yUiIj7H95Z8+fLhiy++wKFDh9CxY0cMA/AQgPhGM0BrAJyRkdfte1etsvbAhScFSCGQQYrV0pqIiPgZp2J8/vnnZnf6FADNAcQ4b6vOKRoAytruP3HiRDOTNVwpQAqBAMmROXNAD0VERCLH119/bZpLLgLAStg9zuvZBnqlWxsAFmyHKwVIQSrWPiRQAZKIiKSSChUqmGkZHD7/7w034A4Aa523seHMLADP2Oa1HT9+3PW9O3bsQOvWrc1s1lDfA6YAKUhdsrVyj1IPJBERSWX169c3O772ORyowlWN9u1dt70K4EsAG1asQK5cuVw74UaPHo1Zs2bhmWeeMV9DmQKkIHXpJBvBx4tWgCQiIoGUOTPwxRfA8OFwODcNMVz6BUAxjihp0AAvvvgi3n33Xde3TJ48GaFMAVKQumILkNJmyxbQYxEREQEDoyFD4PjmG5xyXlXZufzWAsBLL72U4O7ff/89XnnlFaxYsQITJkzAaXvpSAhQgBSkYmOsvQNAGgVIIiISJKJbtTJ1Sdudl3M4R5W8xvcrt/s+//zzqFWrFh5//HEMHWrvsBT8FCAFqThbH6QMOXMG9FhERETsbr77blQDMMN23dMAfgKQWG9t1ieFEgVIQcphC5DS52B8LiIiEhwmT55sltnaAehraypZD8A6AE0R+hQgBakoWx+kTLlzB/RYRERE7Nh9+8yZM/jxxx/RccUKExhZ/ZLyAZjvnCOaPsF3ATNm2HNO8TgTjf2UWLMUTDi0V4JQ9Llz5uslAFmUQRIRkSCTJUsW3HXXXWZEyUpnwfZUZwduGgCgMYAHAGx2XteuXTuz1HbgwAHcdttt2LJlC7Zu3Ypp06aZ2//9918ULVoUwUABUpBKc+GC+co8UlLThkVERAIpp7NOlt377gHMwNuRrJ8FUAnAGgCDAFgNAPr375/oY7E5ZZcuXRAMtMQWpNJcvOgKkG5QHyQREQlS6dKlQ8OGDc35AgULYttdd6FfrVr403k7p4m+A+AHAAWv8ViXLnHdJDgoQApS6Z2/JCzVVgZJRESC2Zw5c7By5Urs3r3b1CW9t3w5Zg4ejHG2+zC7xKDpoSQeRwGSXFP6y/F7ApRBEhGRYJcxY0bUqFEDadL81wkpa548Zofb3QAOOK/LDuAjBlQACnt4HM5169evH5o1a5ZgxlsgqAYpGF26hLTOIX8MkIoqgyQiIiGmfPny5uuPAMqxDxKArs7bWMj9B4CBAD50ayxp4Yy3CxcuIEMGVjOlPmWQgpGtB5KW2EREJBTVqlXLdT5tnjxmaY3LbPuc13FGxCRnc8lSiTwGl+sCRQFSMLL1QDofFRWw6FlERMRbWbNmNTPYGCh9+eWXppB7jjObZB9jy/LujQCGcanO7THOnz+PQFGAFOQB0qX07m22REREQkP37t2xfPlyNGrUCDNnzsT06dOxets2dANwF4C/nffjO93zzmW3O23fHxsbG6AjV4AU9EtsVxQgiYhIGMiRI4dpFGk1gpzHOiUAnxUrZpoiU3Hn9WwbWQjARWfLm0BQgBTkGaQrGd0TjiIiIqErg61shAtoD+7ahUPz5mGJ7T73A9jK2qSZM4G4uIAcpwKkIHQlJsZ1PlYBkoiIhLnCTZuib4UKYA/tw87rsnAn286dQHRgQhUFSEHowtGjrvOOLPwVERERCR/Dhw83Rdxjx441l6OiorBg4UI0mTLF7Ghjg0kWm0yrUSNgx6gAKQhdOMaJNk4KkEREJMwMGTIEMTEx6NOnj+u6vHnz4sEHH8RJwDSYvAlAmXvYGCAwFCAFoUsnTrjOR2XNGtBjERER8YdoD0tn7MTNrFKVKlUwbcECtG/fHoGiTtpBHiClUZNIERGJIH369EmQWQoUZZCCvEg7TTb2GhUREZHUpAApCMWeOuU6ny47R/uJiIhIalKAFIQcp0+7zqfPkSOgxyIiIhKJFCAFIYetk3aGXLkCeiwiIiKRSAFSEIqyddLOqABJREQk1SlACkLRtunFWfLmDeixiIiIRCIFSEEozYULrvOZc+cO6LGIiIhEIgVIQSidc3rxOQA3qkhbREQk1SlACkLpLl82X1mJdMMNNwT6cERERCKOAqQglF4BkoiISEApQApCmWJjzdezUVHIkCFDoA9HREQk4ihACjYOBzLGxZmzF9KkCfTRiIiIRCQFSMHm/HnX/5SLaTVLWEREJBAUIAUbWxftS+nTB/RQREREIpUCpCBz6cQJ1/krCpBEREQCQgFSkDl35IjrfGzGjAE9FhERkUilACnInDt82HU+LnPmgB6LiIhIpFKAFGQuHDvmOu9QgCQiIhIQCpCCOEBC1qyBPBQREZGIpQApyFw8ftx1PloBkoiISEAoQAoyV06edJ2PvvHGgB6LiIhIpFKAFGSuxMS4zqfNli2gxyIiIhKpFCAFmdjTp13n0+XIEdBjERERiVQKkIKM49Qp1/n0CpBEREQCQgFSsDl71nU2Y65cAT0UERGRSKUAKchEnTvnOq8ASUREJDAUIAWZNOfPu85nyZs3oMciIiISqRQgBZk0Fy64zmdWgCQiIhIQCpCCTNpLl1zns+bLF9BjERERiVQKkIJMemeAFMfzahQpIiISEAqQgkyGK1fM17NRUQBPIiIikuoUIAWZjLGx5uv5aP2vERERCRS9CweZTHFcXAMupEkT6EMRERGJWBEbIL3zzjsoVqwYMmbMiBo1amDVqlWBPiRcvHgRWa3z6dIF+GhEREQiV0QGSNOmTcOAAQMwdOhQrF27FpUqVUKzZs1w+PDhgB7XqRMnkNF5/nL69AE9FhERkUgWkQHSqFGj0L17dzz88MMoW7Ys3n//fWTOnBmTJ08O6HGdOXTIdf5KhgwBPRYREZFIFnEB0qVLl/D777+jSZMmruuio6PN5RUrViS5/HXq1KkEJ1+7cvKk63xcRiuXJCIiIqkt4gKko0ePIjY2FvncmjDy8sGDBxP9vhEjRiBbtmyuU5EiRXx+bCULFnSdr1Snjs8fX0RERJIn4gIkbw0ZMgQxMTGu0549e3z/Q7is1qYNcOediKpUyfePLyIiIsmSFhEmd+7cSJMmDQ7Z6n2Il/Pnz5/o92XIkMGc/KpoUeDrr/37M0REROSaIi6DlD59elStWhWLFi1yXRcXF2cu16xZM6DHJiIiIsEh4jJIxC3+Xbt2RbVq1VC9enWMGTMGZ8+eNbvaRERERCIyQGrfvj2OHDmCF154wRRm33bbbZg7d+5VhdsiIiISmaIcDocj0AcRirjNn7vZWLB94403BvpwRERExIfv3xFXgyQiIiJyLQqQRERERNwoQBIRERFxowBJRERExI0CJBERERE3CpBERERE3ChAEhEREXGjAElERETEjQIkERERETcROWrEF6wG5OzIKSIiIqHBet++1iARBUheOn36tPlapEiRQB+KiIiIePE+zpEjidEsNi/FxcVh//79uOGGGxAVFeXTyJZB1549ezTjLQLp/79nt99+O1avXo1gFYjj89fP9NXjXu/jePP9Kf2e5Nxff5Ph93fJsIfBUcGCBREdnXilkTJIXuI/auHChf32+PxD1B9j5NL//4TSpEkT1P8egTg+f/1MXz3u9T6ON9+f0u9Jyf31Nxlef5dJZY4sKtIWkaDXs2dPBLNAHJ+/fqavHvd6H8eb70/p9wT771Ww6xnk/37Xe3xaYgsyTOcyso2JiQnqyFz8Q///RYKL/iYjlzJIQSZDhgwYOnSo+SqRR///RYKL/iYjlzJIIiIiIm6UQRIRERFxowBJRERExI0CJBERERE3CpBERERE3ChACiLvvPMOihUrhowZM6JGjRpYtWpVoA9JREQkIilAChLTpk3DgAEDzHbStWvXolKlSmjWrBkOHz4c6EOTINGmTRvkyJED9913X6APRUQAM36kQYMGKFu2LCpWrIivvvoq0IckPqRt/kGCGSPOjRk/frxr1hvn//Tu3RuDBw8O9OFJEFiyZImZHzRlyhTMmDEj0IcjEvEOHDiAQ4cO4bbbbsPBgwdRtWpVbNu2DVmyZAn0oYkPKIMUBC5duoTff/8dTZo0STDrjZdXrFgR0GOT4MFPqhyOLCLBoUCBAiY4ovz58yN37tw4fvx4oA9LfEQBUhA4evQoYmNjkS9fvgTX8zI/lUjoW7p0KVq0aGGmR0dFReHbb7+96j6qQRMJ3b9Lfsjl6zgz/xIeFCCJpIKzZ8+aujK+2HqiGjSR0P27ZNaoS5cumDBhQioduaQGBUhBgGnZNGnSmLVsO15m2lZCX/PmzfHKK6+YQmtPRo0ahe7du+Phhx82BZ/vv/8+MmfOjMmTJ6f6sYpECl/8XV68eBGtW7c2taK1atVKxaMXf1OAFATSp09vivsWLVrkuo5F2rxcs2bNgB6b+J9q0ERC8++Se5weeughNGrUCJ07dw7g0Yo/KEAKEkzjTpw40exQ+uuvv9CjRw+T/uUnFwlvya1B4wtzu3btMGfOHBQuXFjBk0iA/y6XL19uluFYu8RibZ42bdoUoCMWX0vr80cUr7Rv3x5HjhzBCy+8YP74+Ic2d+7cq/44JXItXLgw0IcgIjZ16tQx2X4JTwqQgkivXr3MSSKLatBEgo/+LkVLbCIBpho0keCjv0tRBkkkFZw5cwY7duxwXd61axfWr1+PnDlzomjRoqYGrWvXrqhWrRqqV6+OMWPGqAZNxM/0dylJ0agRkVQaE9KwYcOrrueL78cff2zOc8zMG2+84apBGzdunGlMJyL+ob9LSYoCJBERERE3qkESERERcaMASURERMSNAiQRERERNwqQRERERNwoQBIRERFxowBJRERExI0CJBERERE3CpBERERE3ChAEhEREXGjAElERETEjQIkEfGbBg0aoF+/fgh1of48Uvv4Bw0ahNatW6fazxPxBwVIIiLiU+vXrzeDXUVCmQIkEQlqly5dCvQhSAr/v2zYsEEBkoQ8BUgiYWju3LmoU6cOsmfPjly5cuHee+/Fzp07Eyy59OnTB//73/+QM2dO5M+fHy+++GKCxzh9+jQefPBBZMmSBQUKFMDo0aOvWqopVqwYxowZk+D7+Mbo/ljJPS7r2Hr16mV+Tu7cudGsWTOPj8X79e7d29wvR44cyJcvHyZOnIizZ8/i4Ycfxg033IASJUrgxx9/dH3PxYsXzfPOmzcvMmbMaI5l9erVCR6X39+lSxdkzZrVPO+33nrrqp8dFxeHESNG4Oabb0amTJlQqVIlzJgxI5H/G8n7N0/Ov6U3z9ly5coV8++aLVs28+/6/PPPw+FwJPs5Jff/y969e3H06FHz/SKhTAGSSBjiG+aAAQOwZs0aLFq0CNHR0WjTpo15E7RMmTLFBD+//fYbRo4ciWHDhmHBggWu2/n9y5cvx3fffWeuX7ZsGdauXev347KOLX369Obnv//++4k+Hu/HN+tVq1aZwKFHjx5o164datWqZY71zjvvROfOnXHu3DlzfwYnM2fONN/H2xlM8I3++PHjrsd86qmn8PPPP2PWrFmYP38+lixZctXzZiDxySefmGP7888/0b9/f3Tq1Ml8X1Ku9W+eHCl9zvbvS5s2rfm+sWPHYtSoUZg0aVKKnlNy/r9weY1BGAMtkZDmEJGwd+TIEaYKHJs2bTKX69ev76hTp06C+9x+++2Op59+2pw/deqUI126dI6vvvrKdfvJkycdmTNndvTt29d13U033eQYPXp0gsepVKmSY+jQoa6fY7//tY7L+p7KlStf8zm5P4crV644smTJ4ujcubPrugMHDpjHX7FihePMmTPmOX322Weu2y9duuQoWLCgY+TIkeby6dOnHenTp3dMnz7ddZ9jx445MmXK5HoeFy5cMP8Ov/76a4Lj6datm6Njx47JPl73f/Nr/Vt685zt31emTBlHXFyc6zr+XF6X3OeU3P8vL7/8sqNevXrXvJ9IsFMGSSQMbd++HR07dsQtt9yCG2+80Szf0O7du133qVixYoLv4XLS4cOHzfm///4bly9fRvXq1V23MytQqlQpvx8XVa1aNVmPZ38OadKkMct2FSpUcF3HJSji8+JSHp9T7dq1XbenS5fOPMe//vrLXOZ9WFtTo0YN1324HGZ/3jt27DDZmaZNm5plOOvE7Iv7cmFSx+v+b55cKXnOdnfccQeioqJcl2vWrGn+f8TGxib7OSXn/wszSFpek3CQNtAHICK+16JFC9x0002mPqVgwYJmCat8+fIJCmsZHNjxzdN9qetauERmr2MhBiHXc1zEZajk8PQc7NdZAUFKn1dSzpw5Y77Onj0bhQoVSnBbhgwZUny81rEl99/SH885uc8pOf9fGCDdfffdyf7ZIsFKGSSRMHPs2DFs3boVzz33HBo3bowyZcrgxIkTKXoMZnj4pmsvYI6JicG2bdsS3C9Pnjw4cOCA6/KpU6ewa9cuvx3X9ShevLirfsYegPA5li1b1nUfPm/WCFl4jPbnzfsyaGDWizVM9lORIkW8Pr6U/Ft6w/6caOXKlShZsqTJQvnqObGwn9lH7WCTcKAMkkiY4e4mLrtMmDDBLOHwTW/w4MEpegzuhurataspWOYSE3d9DR061GQ57Ms0jRo1wscff2wyQ9yZ9sILL5g3XH8d1/Vg9oMFzdZzKlq0qCmU5tJSt27dzH24rMTzvA+Plc/72WefNc/b/m/DRogsYmaWhjvhGDwy8OKyIf/dvJGSf0tv8N+bBfKPP/64KeZ+++23XTv0fPWcuL2fx1yuXDmfHbdIoChAEgkzfDP/8ssvzZZyLl+xfmbcuHFmm3ZKcJfTE088Ybbi802SO8D27NljtsdbhgwZYrIcvA9rlF5++eVEsx6+Oq7r8dprr5kAgLu8mO2oVq0a5s2bZ4I3yxtvvGGWnBioMHAYOHCgCRbs+DyZ8eHOL2ZMGNBUqVIFzzzzjNfHlpJ/S2+wdcH58+dNzRWDmL59++Kxxx7z6XPi8lrp0qWvudQoEgqiWKkd6IMQkeDHLfqsT2HWwcq4iIiEK2WQRMSjdevWYcuWLSbjwAwKe/ZQq1atAn1oIiJ+pwBJRBL15ptvmsJqFjdzizebRbJJoYhIuNMSm4iIiIgbbfMXERERcaMASURERMSNAiQRERERNwqQRERERNwoQBIRERFxowBJRERExI0CJBERERE3CpBERERE3ChAEhEREXGjAElERETEjQIkEREREST0/0CjROjJCwmjAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# get the angular power spectra of the lensing maps\n", "sim_cls = hp.anafast(\n", " [kappa_bar, gamm1_bar, gamm2_bar],\n", " pol=True,\n", " lmax=lmax,\n", " use_pixel_weights=True,\n", ")\n", "\n", "# get the expected cls from CAMB\n", "pars.min_l = 1\n", "pars.set_for_lmax(lmax)\n", "pars.SourceWindows = [\n", " camb.sources.SplinedSourceWindow(z=z, W=dndz, source_type=\"lensing\"),\n", "]\n", "theory_cls = camb.get_results(pars).get_source_cls_dict(lmax=lmax, raw_cl=True)\n", "\n", "# get the HEALPix pixel window function, since the lensing fields have it\n", "pw = hp.pixwin(nside, lmax=lmax)\n", "\n", "# plot the realised and expected cls\n", "ell = np.arange(lmax + 1)\n", "plt.plot(ell, sim_cls[0], \"-k\", lw=2, label=\"simulation\")\n", "plt.plot(ell, theory_cls[\"W1xW1\"] * pw**2, \"-r\", lw=2, label=\"expectation\")\n", "plt.xscale(\"symlog\", linthresh=10, linscale=0.5, subs=[2, 3, 4, 5, 6, 7, 8, 9])\n", "plt.yscale(\"symlog\", linthresh=1e-9, linscale=0.5, subs=[2, 3, 4, 5, 6, 7, 8, 9])\n", "plt.xlabel(r\"angular mode number $l$\")\n", "plt.ylabel(r\"angular power spectrum $C_l^{\\kappa\\kappa}$\")\n", "plt.legend(frameon=False)\n", "plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.2" } }, "nbformat": 4, "nbformat_minor": 4 }