{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Cosmic shear\n", "\n", "This example simulates a galaxy catalogue with shears affected by weak lensing, combining the [galaxies](../1-basic/density.ipynb) and [lensing](../1-basic/lensing.ipynb) examples with models for the intrinsic galaxy ellipticity and the resulting shear." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup\n", "\n", "The setup of galaxies and weak lensing fields is the same as in the basic examples.\n", "The setup for angular matter power spectra matches the definition from the [Matter shell definition](../1-basic/shells.ipynb) example." ] }, { "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.compat.camb import Cosmology\n", "\n", "# almost all GLASS functionality is available from the `glass` namespace\n", "import glass\n", "import glass.ext.camb\n", "\n", "# how many arcmin2 over the entire sphere\n", "ARCMIN2_SPHERE = 60**6 // 100 / np.pi\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", "results = camb.get_background(pars)\n", "\n", "# get the cosmology from CAMB\n", "cosmo = Cosmology(results)\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 window function for shells\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)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Matter" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# set up lognormal matter fields for simulation\n", "fields = glass.lognormal_fields(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)\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": [ "## Galaxies" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# standard deviation in each component of galaxy ellipticity\n", "# this is very small so that the galaxy density can be small, too\n", "sigma_e = 0.01\n", "\n", "# galaxy number density per arcmin2, over all shells\n", "n_arcmin2 = 0.01\n", "\n", "# localised redshift distribution with the given density\n", "z = np.arange(0.0, 2.0, 0.01)\n", "dndz = np.exp(-((z - 0.5) ** 2) / (0.1) ** 2)\n", "dndz *= n_arcmin2 / np.trapezoid(dndz, z)\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", "Simulate the galaxies with shears.\n", "In each iteration, get the shears and map them to a HEALPix map for later analysis.\n", "\n", "In addition, generate the galaxy ellipticities, drawn from the intrinsic normal distribution.\n", "The standard deviation is much too small to be realistic, but enables the example to get away with fewer total galaxies.\n", "\n", "Finally, apply the reduced shear from the lensing maps to the galaxy ellipticities, producing the galaxy shears." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# number of HEALPix pixels in the maps\n", "npix = 12 * nside**2\n", "\n", "# map for galaxy numbers\n", "num = np.zeros(npix)\n", "\n", "# map for sum of shears\n", "she = np.zeros(npix, dtype=complex)\n", "\n", "# simulate the matter fields in the main loop\n", "for i, delta_i in enumerate(matter):\n", " # compute the lensing maps for this shell\n", " convergence.add_window(delta_i, shells[i])\n", " kappa_i = convergence.kappa\n", " gamm1_i, gamm2_i = glass.shear_from_convergence(kappa_i)\n", "\n", " # generate galaxy positions uniformly over the sphere\n", " for gal_lon, gal_lat, gal_count in glass.uniform_positions(ngal[i], rng=rng):\n", " # generate galaxy ellipticities from the chosen distribution\n", " gal_eps = glass.ellipticity_intnorm(gal_count, sigma_e, rng=rng)\n", "\n", " # apply the shear fields to the ellipticities\n", " gal_she = glass.galaxy_shear(\n", " gal_lon,\n", " gal_lat,\n", " gal_eps,\n", " kappa_i,\n", " gamm1_i,\n", " gamm2_i,\n", " )\n", "\n", " # map the galaxy shears to a HEALPix map; this is opaque but works\n", " gal_pix = hp.ang2pix(nside, gal_lon, gal_lat, lonlat=True)\n", " s = np.argsort(gal_pix)\n", " pix, start, count = np.unique(gal_pix[s], return_index=True, return_counts=True)\n", " num[pix] += count\n", " she[pix] += list(map(np.sum, np.split(gal_she[s], start[1:])))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Analysis\n", "\n", "Compute the angular power spectrum of the observed galaxy shears.\n", "To compare with the expectation, take into account the expected noise level due to shape noise, and the expected mixing matrix for a uniform distribution of points." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAioJJREFUeJztnQd4U3Ubxd9udpmyp7IVkClLUBBEliCKfMoSUVBkCgKCCKgoKCAKIiAgioioIA5wsJcge8oQEGRvKLO0+Z7zb+7lJk3apEma5Ob8nidwc5Pc/NM2ybnvOG+YxWKxCCGEEEIICXrC/b0AQgghhBDiHSjsCCGEEEJMAoUdIYQQQohJoLAjhBBCCDEJFHaEEEIIISaBwo4QQgghxCRQ2BFCCCGEmAQKO0IIIYQQkxDp7wUEO4mJiXL8+HHJmjWrhIWF+Xs5hBBCCDEZmCVx5coVKVCggISHpxyTo7DzEIi6woUL+3sZhBBCCDE5R48elUKFCqV4Hwo7D0GkTvthZ8uWzd/LIYQQQojJuHz5sgoiaZojJSjsPERLv0LUUdgRQgghxFe4UvLF5glCCCGEEJNAYUcIIYQQYhIo7AghhBBCTAKFHSGEEEKISaCwI4QQQggxCRR2hBBCCCEmgcKOEEIIIcQkUNgRQgghhJgECjtCCCGEEJNAYUcIIYQQYhIo7AghhBDiFTp16iSPP/64z5/nzTfflEqVKgXMcQIJCjtCCCGEeIUPP/xQZs6cKYE6Z3XBggU2+1599VVZsmSJmIlIfy+AEEIIIeYgNjbW30twiyxZsqiLmWDEjniFGzduyNatWyUhIcHfSyGEEOJjvv32W7nvvvskY8aMkitXLmnYsKFcvXo1WSq2fv368sorr0jv3r0lR44ckjdvXpk6daq6b+fOnSVr1qxyzz33yKJFi/THIOKXPXt2m+dDpA0RN2f89ddf8sgjj0ju3LmVuKxXr55s3rxZv71YsWLq/1atWqnjaNftU7GJiYkyYsQIKVSokMTExKjbFi9erN9++PBh9fjvv/9eHnroIcmUKZNUrFhR1q1bJ4EChR3xmLNnz6o/7Pvvv19Gjhzp7+UQQgjxISdOnJB27drJc889J3v27JHly5dL69atxWKxOLz/559/rgTXhg0blMjr3r27PPnkk1KrVi0lvho1aiTt27eXa9eupXlNV65ckY4dO8rq1avlzz//lJIlS8pjjz2m9mvCD8yYMUOtX7vuKJX8wQcfyPvvvy/bt2+Xxo0bS4sWLWT//v0293v99ddVGhcBjVKlSqmfx+3btyUQYCqWeER8fLw89dRTsm/fPnX9m2++UWdAhBBC3Kdq1apy8uTJdH/efPnyycaNG126L4QRRAzEXNGiRdU+RO+cgRP/IUOGqO1BgwbJu+++q4Re165d1b433nhDPvnkEyWkHnjggTSt/+GHH7a5PmXKFBX1W7FihTRr1kzy5Mmj9mMfXqszIOhee+01efrpp9X19957T5YtWybjx4+XiRMn6veDqGvatKnaHj58uJQvX14OHDggZcqUEX9DYUc8An/c+KPXgMC7efOmCmEHIkeOHJHo6OgU39je5NatW+r5CCHEFSDqjh07JoEMhFqDBg2UmENECxG3Nm3aqFSrIypUqKBvR0REqNStUQgiPQtOnz6d5jWdOnVKicfly5er46AsCBFAfOa7yuXLl+X48eNSu3Ztm/24vm3bNqevKX/+/Pr6KexIUIOQ9oQJE2z24c30999/qzd+oIEQPeo9kC5AGgBnZUgfexscHx8uOPND3QjO/L766iuvPw8hxHyk10mnJ88Lcfb777/L2rVr5bfffpOPPvpIpSbXr1/v8P5RUVE211GjZtyn1c6hvg2Eh4cnS+siO5QSSMOeO3dOpVIRRURwoWbNmurk2hektH5/Q2FH0gTewN26ddOvV6lSRTZt2qS2d+zYEZDC7qefftKbO+bOnasuONMcOHCgEnwpFea6AlITKCiGoNN+FmDOnDkyePBguffeez1+DYQQc+NqOtTf4PMSkSxckEqFmJo/f75Xjo20KWrj0GCROXNmtQ+1bCmxZs0amTRpkqqrA0ePHlX13/ZiLKUGv2zZskmBAgXUsdB8YTx29erVJVhg8wRxG9RXoLZCOxNCISy6iDR27twpgQjqH+zB2SZqM2rUqKG6nNJyxhUXF6fOEtHZhQJao6gzijtCCDHLif0777yjRChSnfjsPHPmjJQtW9Yrx8fnMbpNcUL8zz//qIxHat54aJb44osvVDMH1vfMM8+ojl0j6ISFZx3S3RcuXHB4nP79+6u6Opz47927V534Q1T26tVLggUKO+IWqJ+DqEMdAqhbt64qKjXWSyBiF8jCLjIyUj7++GMpXry4fhs6pJ544gkpV66cfPbZZ+p1uiJw8cFTuHBh1cr/77//6rchxTt58mSVUtCEnbOOMUIICSYQ2Vq5cqWKjqEjFLVt6CRt0qSJV46fM2dO+fLLL+WXX35R3y34/EytKQ+f2xBrlStXVh22PXv2lLvuusvmPlgjUsj4zHZWhoPH9e3bV/r166eeG1YnCxcuVMIxaLAQj7h06RK+rdX/ZicxMdHSpUsX9XpxKVSokOXkyZP6bbGxsWp/kSJFLIEG1pc5c2a1vpIlS6p98fHxljlz5lgqVqyovybtUqBAAcuYMWMsly9fTnasXbt2WZ577jlLdHR0ssc1adLEsmTJEvV8oGHDhvpt69atS/fXTQghJLS0BiN2xGXQjo6zIpAhQwZlGKl1M6HeQqshQ2j+0qVLEkigYwr1GgApUy1yh8aGLVu2qCYH1NlpICKJkHyRIkVUUTAej4YItM2jrX369Ol6Khp1GyjcRaQSZ5hI7Wr1ev/73//0YzIdSwghxNdQ2BGXQNjdWGMwbdo01TBhxNgcsGvXLgnU+jpN2GlAhD366KPKtgXGlnAm17h48aKqJUFBLVzGf/75Z5t0xIABA+TQoUOq/sNRcwTS1pr1C2o2AsXAkhBCiDmhsCOpgggcPIo0UYLaAxSm2mOsswu0BgoU4DoTdka0Jordu3ercTdaS7uxqQL1GajVQNcVimwLFizo9HgYbaN1aSHqZ/T8I4QQQrwNhR1JERg8Yu4fOp4AZvHBNdwRxoiVTxooIK7gyI6u04ULEQIT+eUXhBNFtmxBWA7qCYtG8ajLETtHoLsL6daDBw+qQloU8yJCiYJeiETsQ8TOFZiOJYQQkl6EodDO34sIZuBUjagMaspc/aIPFvCngcicJkZKlCihukchchwBc0iMiQGoV3MrOnXjhgjc1nH5778728Z9J07ALM6146EbNUsWkaxZ1eXA6dNy5Px5wdTAek8+KdkrVxa5+268qKT/7QZOe5Pr16+rWkT4MuFvBJE71CgSQggh3tYaFHYeYmZhN2bMGFVDBmASifqz1Ex2UYsGGxCMjEGUz6npL/7s/vhD5OOPMRJC5Px58SsYhWMUesb/CxWC1bpHh0dzxaxZs9Q2Ur3GOr70BrYs48aNU00eGG5NCCEksKGwSwcwDBgXuFhjPqrZhN2vv/6qasO02jJXxQjmBsL0F0DgJRtTg87UL74QwSiyPXtcXxAGOKOWTbtAbGXKBHdgkStX7vyvXQzXLdZLmusOUGdXrBjysyKVKiVdMFkDPnguTqvAzxMNGgD1ivPmzRN/0bJlS+XLBOA/BXNpT6duEEII8R0UdumIGSN2qEerVq2a6ggFGBczfPhwlx6LxoqxY8eqbQg81OQpYN47caLI1KloNbV9EEwkYf4IsWYUb5qAw4Bla2dpWtBSxJlEpFndujIX60MzxcGDSRdt++jRpDo+V8HvGwIPF03wlS8PL5hkd0XjCaKZiGIiDYt0rD/+Xm7cuKEGdeN/DaTbYWOjde8SQggJXq3BWbHEBtSBIaKjiTpsDxs2zOXHG1O1O3fskEcgFj78UGTBguSiqW5dEViotGwJUznxdUfsNTiaQ3hVrZp0sQe+dBCgmtCzF39WHzydy5dFVq1KumggZVumzJ2oHv6vUkUic+aUp556SkV5IargAdihQwdJb1avXm0j6sDs2bPlv//+U3MeIfoIIYTY0qlTJ/W9iM/uQIfCjiSLuMHqQ+sMRV2YNhbLFWB5grjP0/Bwe/vt5LVz0dEi7dolCTonI128jcsdsVgbIoeORsdAlB46hEnUItu2Jf2PC6J8RjBgGh5+uMyefWd/6dIytEQJQevHnyLy9Zdf+kXYaWly7YMK3npo7lixYoXUqlVLGSwbR60RQkgwsnz5cuU9ijFj2d1ojjt8+LD6DIRxfSWcmFvBPPBgSXBS2BEbNANeDE/+4Ycf3EsXHj8uFb/9Vo4gu4rrRlGHWrvu3UVefFHEOq0ivTAKu7vREJEWIG7xWFyeeOLOfrxGTehp/0MYx8fbPn7vXsm7d69Mtl69+vvvcqt2bYlG1LJGDZEHHkhKOaeTsENNHZpjunfvLs2bN5fTp0/L33//LQ888ID8+OOPUr16dZ+vhRBCgoXY2FgJGnw+4MzkmG1WrDbvtXTp0u49cNgwiyUyEuczNpfEatUsli+/tFhu3rT4i/bt2+vzWnfs2OH7J8Rr3brVYpk502Lp2dNiqVHDYomOTvazSXbBjN2nnrJYxo61WNautVhu3PDqsk6cOKH/HKpWrarvP3jwoKVMmTL6bRkzZrR8//33Xn1uQoi5SEhIsLzzzjuWYsWKWTJkyGCpUKGCZd68eWpOdoMGDSyNGjXSZ2afO3fOUrBgQcvQoUPV9WXLlqnPmp9++sly3333WWJiYiw1atRI9vm8atUqS506ddTxMZv8lVdescTFxem337hxwzJgwAB1G2Z333333ZZp06ZZDh06lGyOd8eOHdVjFi1aZKldu7b6rsuZM6eladOmlgMHDujHtH9cvXr11H48vmXLljbPjfXkyZNHrR/H3LBhg3679hr/+OMPS5UqVdTnas2aNS1///23z7UGhZ2HmE3YRUVFqddz//33u/6guXNtBMrtsDDLHBHLAyKWA/v3W/wN3kzam/Tq1av+WQRE2p9/Wv4bMMDylYjlYGoiD5cMGSwWfKgMGWKx/PqrxXL5skdL+OKLL/Sfw+DBg21uO3/+vPoA024PCwuzjBs3zsMXTQgxK2+99ZY6IVy8eLHln3/+scyYMUMJnOXLl1v+++8/S44cOSzjx49X933yySct1atXt8THx9uInrJly1p+++03y/bt2y3NmjVTIvHWrVvqPhBbmTNnVp9D+/bts6xZs0Z9L3Xq1Elfw1NPPWUpXLiwOhHFGiCivv76a8vt27ct3333nXqOvXv3qpPaixcvqsd8++236rb9+/dbtmzZYmnevLkSlxCqAOJME2R4HESpI2HXs2dPS4ECBSy//PKLZdeuXep2vGbt/tprhGDFzwT3qVu3rqVWrVpp+nlT2KUjZhJ2eENpX+w4S3KJ//6zWHLkuCNGevWyjOnVSz/O/PnzLf4GZ1RYC84Y/Q3OYMuXL6/Wk1fEcvLTTy2WgQMtlvr1LZZMmVIWehERFgsibX36WCz4uZ45k+bIJT5o7MEZ6LPPPmtztoozUnxIEkKI8bMiU6ZMlrXILBjo0qWLpV27dmr7m2++UZG2gQMHKoEGcaahiR6IMA0IIkS15iJQYD3WCy+8kCyCFx4ebrl+/boSbDjG77//7nCNy6zPceHChRRfy5kzZ2yyOVq0D6LPiFHYIWqIIMjs2bNtvj8h9EaPHm3z/BCIGj///LPah/X7Umuwxo7YjA/TyASPuNRAQ0GnTiIXLiRdb9tWZNw4KQKPNnTCWmfGYiSZP1vEtXForowS8zWobWvXrp3yjzslIjPOn5eBo0Yl3YipGmi6+PNPkbVrk7pt0bBhbMzYuDHpMm5c0j546z34YFKHMS5Fijh8XpzEafV1MJuuWbNmsvvA7gTNMigcHjlypNr30UcfSXR0tLz//vve/2EQQpKDjn2MTkxvUAeNzxYX65bxfaHbWVm5deuW3G9tinvyySdVpz1GUH7yySdS0kFTmvFzCBONSpcuLXus/qbbtm2T7du3q6594+cYvFUPHTqkxlZGRERIvXr13HqZ+/fvVxZe69evl7Nnz+perZiJnpoBv9FpIT4+XmrXrq3vw1xx1CZr69eoUKGCvp3fWkeNmuYiTj6rvQGFHUm7sIPJMKZHAHjOffKJMuxFZ6xPZ8amweokUIQd0IQdwLi2gQMHJt0AyxfNFw9NJgCj1DRLFczEhfAzgg8RXD79NOk6jJQhpJ96Kqkhw2o8jN8DvPMAOsUg1pwJTxgWQ9x17dpVGXB/8MEH6neK6RmEEB8DUYcxigFMHAzgrc12BfHZb0Dzw8T3yaZNm5T4gphKy3O8+OKL0rNnz2S3QRQZm+LcoXnz5lK0aFGZOnWq8haFsIOggyj1BRB8GpoRvCYmfQWFHdG5avBpS1XY7dwpogkS8PnnSWO5rAIKwgFvFETs/IlXOmK9DGbu1qhRQ50x4ox0165dUh7+eo6AQTPsYXAB586JrFmTJPIg9jZtSorkaRw+LDJ+fNIFZ4QQeG3byq9Ll+p3adSoUapr7Ny5s/K7e+mll9T1F154QZ1xwxKFEOJD7Kf1BODzlitXTgk4RLmcRcxgnQWrrEWLFqkpRk2bNlVjDI1gTKUWuYItCaY4wWYLVK5cWVlvOTshx8kmBBKsmho2bJjs9mjryStOTo1m9Xv37lWiri4yHFZvz9QeZw++S3C/NWvWKJEIEMHDLPXevXuL33E70UtMW2O3bds2vbaqa9euKTcCVKx4p/ard+9kd6lYsaI6TkREhKrH8Bfo2tJeE2o+AoUPP/xQX9frr7+e9gNduWKxoMbkjTcsloceQveLw/q8YxkzWt4RsVQUsfy9Z4/Lh+/evbu+zrx581qOHDmS9rUSQkwDPrdy5cplmTlzpmp02LRpk2XChAnqOrpd0aWKfWDQoEGqcxVNWsb6M9QbowYN9W0tWrSwFClSxHLT6qCA7yPU3L388suq3g01egsWLFDXNdBIgeYJ1HKjux/HnWut0UMDB5rAsJ7Tp09brly5ohoksGbUEqN5YsmSJZZq1arZ1IOjwQPPi+aQkydP6k0X9s0TvXr1UjV16LI1Nk/Yv0ZjjR9eB/ahjs9d2DyRjphJ2K1bt07/EscfrVMGDLgjGsqXt1gcFII+88wz+rHwBvUXzz33nL6OzZs3WwIFdFuhCBjrKlGihG4L4DH4EJkxw2J59FGH9jPKgqZUqaROWxcEHgqCH3roIf1niK40o90AISQ0wWcWul5hjYVGAjSpNW7cWDVm4SQQJ9XGzxFYfqCL1Sh6fvzxRyXuIALRNWv/XYEO1UceecSSJUsW1YABS5W3335bvx1NCH369LHkz59fHeOee+6xTJ8+Xb99xIgRlnz58imBp9mdoNkC3bjo4MXxsF77Rr+pU6cqwYjPaGd2J3huNJflzp07RbsTCrsgxEzCDmcv2hc4zrAcgm7KsLAkkQBvNvi1OeDdd9/Vj/UlfOz8xIMPPqivI9B+Rw0bNtTX9ueff3r/Cc6exSeU5cz991tuO+u0Rev9Z58lRf6cHuasEp/aWmFd4DUhSggJOVztWCVp0xquz4oiIVVjh87JZFy6JIIxWNpYlbfeSir0T21mrB/r7LQau7vuusu9KRrpwP/+9z99+6uvvvL+E+TKJfL88zLqoYcEvVjd0I2FWj5rAa8C3bdduiTV1+B/XLcbm5MrVy5ZuHChZM2aVV2fN2+e3jVLCCEksKCwI653xfbogZ7wpG0UzPbt6/RYRmHnr85YvJ7jx48HVEeskVatWumFupjZeht2Jz4ANicwfJkaHi6RaLrAzwR2KcbWfoj66dNF0L5frpzImDEi1i5agOYOiE+tq2vYsGHy3Xff+WS9hBBC0g6FHXFN2H3zjciXXyZtI/KFLtiICKfHQqeTFuHxV8Tu4MGD+nYgCjsMpkanGIAVCYZWexsIW+3nX61aNeUVpaJz6Nzavl1kw4YkaxVjNPPvv0UGDEjqyG3dWmTJEhXFa9asmYzSPPcEwdsOshWzcQkhxA3q16+vPOnwGUi8D4UdSd3uBJ5K3ZDIszJxooi1xdsZiOxoUbt///1XGQWnN4FodeLI086X6djff//duc0Jom/VqolMnixy4kSSWIfZsQYiiPPni8BKACncSZNkwEsvybPPPqufCLRo0UL3xyOEEOJ/KOyIw4idXmPnaLrEM8+4dDyjUTG82vwp7AIxYgcQBcuSJYva/v7775V3nDfRpk2k6l8HIY/6yRUrRPbtExk0CDbpd26HCfLLL0tYoUIyIzZWnrC6qR89elQ5zCfNziaEEOJvKOxIyqlYJ9MlXMHfDRTBIOwyZsworZHuVL0pl5SZp7eAeacWsUNaHKbILoHRP++8g1BrUgreauSpuHxZIidOlG+3b5clMTHSBAafq1bJ5s2bvbZuQgghaYfCjjgXdilMl3AFf48WCwZhZ5+OxYgxb4FZi9qcXDi+G0fbuATu/+STSVMuUEv3/PNQovrND9+8Kb8gmAdH9/feE7l502trJ4QQkjYo7IjDGrss+FJHLZX2ZY1i+wYN3DqevyN22pzYHDlyJDUNBCgNGjSQPHnyqO0ff/zRa/WILqdhXQG2NlOnJs2uHT06aSatldI4/rx5IsWLJ90GWxxCCCF+gcKOOIzYFfnsM4R8kq6gcN7QDekquXPnlnzW+YOI2KVnHdbNmzfVHMNAj9YBRNJQpwZQY7dgwYLAE3YaEMj9+yMcKonz58uqSMO4aTRgvPZa0oxaRHpxnRBCSLpCYUeSCTv0ReaCpxmAz9rs2SIZMqTpmFrU7uzZs3L69GlJLw4fPqxqzAK5I9aZWbE30rGIvmrDrYsXL+79n0FEhIQ//riMbdZMqosIHO0sWu0lIo5IzSKq98ILIvv3e/e5CSGEOIXCjtiIAbiZzYIThgvTJQK5zi5Y6us0atasqbz/ABoetNq4tLJy5Uq5deuWHq3TjIW9Tb169eQvEWkjIl8OHpxUh2c1XRY8P9K3ZcqItG8vsnevT9ZACCHkDhR2xCZi966I6A51qUyXCOQ6u2ATduHh4XoTRUJCghrb5Qm//vqr99OwDnjQ4Hv3A4yNIeQOH04yONZMjxE5hbk1JlqgbhP3I4QQ4hMo7IiNsGtu3bag+zGV6RKuwIidf2bHavV1ERERqiPWV1SsWFFiY2P1KKGqo4T/HVKxqHF8++2kujxN4CGtD4EHL0R44xFCCPEqFHbERtglDQETCUNaMJXpEq5Qrlw5PQ3IiF3qIhg/L7BmzRo1sSMtwDR4j1U0wbvOl2N7IBzr1KmjtpE+/tsYjYPgQ3oWETw03+TKlbQf4g/CFU05iFIyRUsIIV6Dwo7oXI2Lk8x3Rk945ZiYYFGiRAl9+oTW0JBeVieY6nDXXXdJMAAB7I0mihTHiPk4HbsCkyvswcxgdMkeOiTy7ru2Au/rr5MEHposYKVCCCHEIyjsiM7ta9ck0svCzlhnh+YMdKv6mtu3b8shiAhrR6yvGgd8wdNPP61vv/nmmzJy5Ehl3eJ3m5NUGig0kI51CgQe7FDwN4BUbe7cSfsTEpJq8zDxArV558/7fM2EEGJWKOyITpjBx84Xwi696uzgXwdxF0xpWA0I0ZYtW6ptCLo33nhDKlSoIEuXLnXp8Wi80CJ2qH2rVq2a+JrKlSvrs4URsUvVrxCzcSHgIL5HjrzTZIE5uWPGiCDCi5FmBsNsQgghrkFhRxTqy9j4RepFYWdsoEiPOrtgrK8zMnv2bOnfv7+qXwP79u1T0ynat28vp06dSvGxW7ZskfPWiBceE2k0EPahwXLt2rXV9vHjx/U0eKpA4A0Zgry5SL9+IjExSfsxueL116FyRaZMQQjWh6snhBBzQWFHFPHx8ZLBWP8WxBG7YBd2iH6NHj1aNm/eLLVq1dL3f/nll1KmTBmZPHmy01rF9E7DOqqzSzEd6wikZN9/P8nIGD544daPJYjYF18UqVQJ/i1eXjEhhJgTCjuid8TaSDkvCrtSpUrpA+gZsXMdpGBXrVolU6dOVfNuwcWLF6V79+5K8G3dujVghJ2xzs5hA4UrFC6cVGu3a5fIE0/c2Y/rjz6adPHDzGFCCAkmKOyIz4UdRB0iTWDv3r36RARfYUwFBrOw04yLn3/+efVz69ixo75//fr1UqVKFenbt69cuXJF7cP/a9eu1V83RomlF6jly2AdO5dmYaeBv5VvvxXBSLTqGFhmBVE7TEFBB+3Jkx6umBBCzAmFHdE7Vn0l7Ix1dmhqgEhJj4gdhEaBAgXEDOTJk0dmzpwpy5Yt00Uy0rHjxo2TsmXLynfffSfLly9XKfX0jtaBmJgYNRYNwH8vrR58NqBub926JM8767g1ZXKMqF6pUknp2zSeJOBvcMKECek6v5gQQtIDCjvi84hdeo4Wg9jRInbwz0PEy0zUr19ftm3bJm+//bYeITt27Ji0adNGOnfurN8vvYWdx3V2zsDvDybGMD6GyTEsUwCilP3744xBZPFitw6JzuHGjRtLr1695OWXX/bOOgkhJEAw17ceCVhhl16jxSByNN+3YE/DOiM6OloGDx6sDJ+bNGmi7z937pz6H52wDz30ULqvyyt1ds7AiDuYHCMa27Ur3JyT9u/bJ4KfQYsWSbe5wJ9//qlHFLXUNSGEmAUKO5Iuqdj0itiZpXHCFRCR/Pnnn+Xbb7+1STmjsSKb5g2XjjzwwANKdLoq7GCxk6rnnT2YIgILlI0b8ULv7P/xx6QJFhhhlor/HX5mGrBnwd8+IYSYBQo7ki4Ru6JFi6rxXr6O2IWSsAOYqvHEE0+oGa2vvfaaStWOgcmvH8iYMaNUtzY74PcA0eQMRFVbt26tBCnm4rpN5cpJzRWzZ4toohb1dkjX4iTil19cEnbg4MGD7j8/IYQEKBR2JF2EHQSIFrXDWDGtk9PbmKkj1h2yZs0q7777rmqu0MSVP3C1zu69996TBQsWyMmTJ2Xs2LFpezKkYzFbF804gwYhR520HyPLmjYVefJJhORsHnL06FHZvn27zT6XDZWtzJo1S4oVK6aaLwghJNCgsCPpkoq1r7NDfZgvCLWIXaDhSp0dOlLR/KGxZ88ez54UkWCMIINgM9YWwjIFHcQffZQ0j9ZBtM7+b8YVML8XNXoY9+Z2KpkQQnwMhR1Jl4hdetXZaV/SaCAoDMNbkq6gvk8bheYoYoeu5RdeeMHGy3D//v26TYtHlC4tsmQJQmpJ0ywAIsM9e6IAUGTzZo+FHTpqEXEGly5dUobRhBASSFDYkXQTdr7ujEX0RPuShjlvesxJJbagjrJq1apqe/fu3XLmzBmb22fMmJFM8MHb0N10aIrp2fbtk9KzGE+msXGjWKpVk0aLFgkqPbNnz54mYYe6QaxX48iRI95ZNyGEeAkKO2KaiN2pU6f0DkemYSXg6uzw++kP7zkrdevW1bfR/OFVcuZMMjJetSqpWxaaLzFRXklIkN1Ip1apoo9pc0fY2Rsve8WImRBCvAiFHUm3GjtMT7gLdhU+itixvi7w6uyMwq5Pnz5y4cIFtf3MM8/YmAN7XGfnjDp1VAoWNXi3rCliJOh7LFki8y0WKWSNumneh6lhL+QYsSOEBBoUdiTdInbGdCxSdN4e50RhFxjUqVNHdUEbGygWL14sc+bMUds5c+ZUnbDaaDSfCjsQHS2WgQOlQd68ssiwu97Fi4IWnucsFjnkouUJhR0hJNChsCPJhJ0FkQ3NOsKH6VhvR+1C1eok0IiNjZVKlSqpbViL/Pfff9K9e3f99vfff19FbkuVKqULQK+nYu1Avd/q48flMRF5CycX+fKp/bBxniYiuZ55Bqot1eNQ2BFCAh0KO5IsFWvJlOnOyCYfNlB4u87OGLG7++67vXpskrZ0LBpaWrZsqXeSwkC5U6dOuqExmlw0Yeepdci+ffuc1sv99NNP+nbWLl0QIpT9hhq/PFu2JBkbf/opFu30OVhjRwgJdCjsSPJULISdj/BlxE77Ug8PD1cGsiQw6uw2o8ZNRGJiYuTTTz/Vo3RAS8fCsBpzftMCUvrt27eX0qVLqyggnsMeo81JU5gXZ88up0eNkkdhWqzdEBcn0q2bSMOGIocO+Txih3KE69evp/nxhBDiCAo7ERk3bpyUL19eypUrJz179gxJ01EbYWcd/eUL8HP2RcQOvzP4oYEiRYooEUH8h7HjVWPIkCFKeBkpW7asvu1uOhaeeFOmTFGC7ssvv9T/Drp162YzFeL8+fOydu1atY37aml6/P8rTjZEZFEhtFFYWboUoWWRiRPxJPpuHNte2J04ccLGk89V5s+fr8bswWvx7Nmzbj+eEEKcEfLCDmfNH3/8sWzatElFkPD/n3/+KaEs7MJ8KOzgc6al3yDs8OXsDfDlDcNYwPo6/5MrVy6b6CxOmgYMGJDsfkZh504DBWr30KTx4osv6ibBSO1q9OrVS5+Z++uvvypjYT1aZwV1fpkzZ5bLItIzQwaR337DWUHSjbDN6dFD5OGHUbypdkGA2UfYIPbcjTTiBKRjx47qWOfOnVONJYQQ4i1CXtgBGI7euHFDud/jollyhBI34uJEi3GF+agj1r7ODnV93qpRYkds4NGiRQv1P1KviKxFO2jIcbczFgINXniVK1eWdevW6fuRikUdH8Z8aUBIvvXWW8nTsFawLu1vBY+9jXFkiCIjHauBrt4KFdRYsn+dpGfdScdCzLVp08ZmVrIWaSaEkJAQdvDBat68uRQoUEB9EGNwuD0TJ05UNVUZMmSQGjVqyIYNG9zyVnv11VdV+g7P0bBhw9AsvLca+yp8LOx8UWfHjtjAA6lX2JosWbJEateu7fA+7qZikWJFV60WgUNqdenSpTJr1ix1QjZ8+HCbObRDhw6Vr7/+Wm1ny5ZNRfmMaH8rOLlTAi1rVpFPPkkaTabVaV67psaSFe7aVQpaH2c8+XPn5ASlHog22jd9EEJIyAg7RHUqVqyoxJsj5s6dK3379pVhw4apIm3ct3HjxjYeabBegJiwv2A8EAxT0TGHM3akVFCL42jGpelJR2Hni85YRuwCD6RGYUr8ECJhToCnnSaSXInYLVy4UN+GiNu2bVuy4w8ePFiJPw1NBDZq1ChZ1ND4t2LTUYsULE46DCbKebdvF5yGPGVXQ+hqxA7ic9q0afrPBk0+gBE7QkhICbsmTZqodEqrVq0c3o6IQNeuXaVz586qjmfy5MmSKVMmmT59un6frVu3KgFhf0GE7o8//lAf7viCwYctUjWhWGMXhqiEHyJ2vhB2IRlxDWK0dOzJkyf1ejlHQKD99ddfahtNB0i7OmuS6devn6qdNWJMw6Yq7ABqTXGM338XKZgUq8MQsrnwwjt8WLK7Iez27t1r4+WHzymtcxsRu1Bs2CKEhKiwSwl0o6HZAelTDZwF47qx/iYl8AWBKB1q7PDFsXz5cpXecQZGD12+fNnmEuzgSyXcWBTuY2GHzsioqCivpmKNX8olSpTwyjFJ+uBqA8WuXbv0WcAouUgNjCybOnWqOmGD2G/durV7wk4Dny87dshKq7gDZTZtEiRUH3ZR2GEdaFACzz//vHTo0EHvEEa9HeboEkKIhLqwQ5caxFjevHlt9uM6zv5d4YEHHpDHHntM7r//fqlQoYL6AtCKvh0xatQo5ayvXSAMgx0I5IzGiIGPhR3SYZp4Rl1VWuwi7NG+lAsWLKgitiR4cLXObv369TbvW1eAiMJnASJmqLFLk7ADOXJIrzx5pJ2IJE27TZo5u0RE2qGmNxU/Os2gGbz++uvq/5IlS+r7mI4lhHiLoBZ23gLF1ogUICKA4myjgao9gwYNUrYa2uXoUd3eNGhJrzmxjtKxKFr3tHgcvwfY1gDW1wUfrnbGGoWdKxE7DQi6CIzJcwDKMbR0borCztokgTaMpvC8a9BA39/xwgWxVK0qgukVToDfnUb+/PnV/0ZPPzZQEEK8RVALu9y5c6sPbPs0Bq7ns86C9Db4EsAXhfFipnFi6SXsvNlAwY7Y0EjFasIO73nYnXgDlG5oNZkHDx506quIdCkarUA07v/bbzKpdGm5Yb09bPduqE0ReOc5OAYatQBqeTUhaRR2jNgRQrxFUAs7pPSqVKmi7BQ08MGM6zVr1vTr2oIJf0bsvFFnR2EX3BQqVEgZBaeUikUtKyLqACUT3ky3a38zqJ91ZjZstDTBxAgJD5ctdesK5GXSwDQRiY+HeZ5I48YI0dnUsGoROy1aZ5+KZcSOEBIywi4uLk51teICDh06pLa1gmVYnaAw+fPPP1dn++g8QwQKXbIkcIWdNyN27IgNbhA102ouETVDI5M9Gzdu1DtHXa2vcxVX6uySCTvr6DrEF7GafWjM0Eo4/vgjydT4p5/UVXT6QjRqqV8NPF6zX2HEjhASMsIOH+hobMBFE3LY1hzm27ZtqzyrcB1+dRB9GNFj31BBnAMhnCmdhR2+HLUojacRO3rYmScdi4i7I5GT1vo6Xws7EC8ii+vVSxJ0WkQO81+bN1fGxicNjRPGiB1SytqJCF5zSuP1RowYodLP2sxbQggJWmFXv359daZuf5k5c6Z+nx49eqgPXpwV4wvA2x/8ZscfETtEabR0LKKwiMymFUbszN8ZG8jCDqgMAkyNMVXC2FX/0UdSqE0bKetA2Bnr7PDZ5awRC93/b775pmzZskVGjx6dlpdICAkhAl7YEXMKO/s6O61+Ki1oX8aYYGCGZpZQJKXOWJzIaabhsBgyNh14A+PJgDvCTvtfjF52uXOLYOwhJuVYmySyHjwoG0XkBQg7u6YuVyxPYJWipaH/+++/tLxEQkgIQWGXRjDiDJMuqlWrJsGOv4SdN+rssHat45BpWHN2xkI0aZ3v1atX10dxeQtE3iIjI5M14jgTdpp3JTwTNWskG5Ni7HvpJRFMyShfXu1CqcOnIvLEnDki1u5aVy1PjJE846hEQghxBIVdGoGr/e7du/URR8GMP+xO7CN2v/32W5qOgWJ7DQq74AW/O81rzj4VmxZjYneAqCtevLgesXM03ksTdrBRypAhg9qGbYlmq2QUfjo4cfnrL1mNRgorhWBmDKsW/O+i5YlRNMKvkePHCCEpQWFH/BaxQ61U9uxJEze/+eYbWbp0qdvHYH2dOUB3qPb7w5QIYyOBL+vr7E8KcJJj74uJ+jfNrsSYfjXW2WG6hdb5akPGjPJx2bLyuIic1/ahmaJOHZHx46Wk4WTElYgdOoa1sWqEEOIICjviN2GXJUsWee+99/Tr3bp1c2h1kRLsiDVfOvb69es2ETCtvk5LxfqClBoojMLKXtgZrzurf0OpwA8iUklEEjRhCs+7Pn0k/8svS0GrJ58zYWc/i5bpWEJISlDYEb+lYrVZnrVq1dJTUe+++65bj6ewM3dnbHx8vGzenGQBXKJECcmTJ0+6CztHjRMayTpjHaBF+y7HxkrEqlVJJsZWwn74QTbEx0t1a3c4Xq899t2y2vg8QgjxWNgxBWD+iJ0Fhd8ZM6bbc6MQ/tNPP9WL10eNGqVSca5CYWfuztjt27frUVxf1Nd5U9g5rLMzCDtldRIVJYIoNcyLc+ZU+wvEx8tqEemZkCCHDDWjGozYEUJ8JuwefPBBfRsmwMR8wi4RheGag346gSaKV199VW3funVLpWRdLRDXvoRz5Mih5nASc3XGpkd9nX19pn1nrCcRO8yY1U6IbTzsmjYVwTQda7Q6SkTGiki2jh1FzuvVeCqCpwlDDUbsCCFeE3bGgubBgwfb3NagQQN3DkUCNWLnxRmc7jB06FC9M3H58uXyxRdfpPoYo6kro3XmithpqVhjfZ0vhV2xYsV0GxV3InYOvewMGEWZvTmxwDZl+XLZCZFnJR+ELKbsrFunrmN2rf1JjiNhx05ZQkiahJ3m2eTog+S84SyTBHGNXTrW1xnBUHd4A2r069dPzp07l+JjYNyqnWxQ2AU/MJfWZqnaR+zQNYuRgb4C1iVa9A21nsbPN08idprHokNhB6KiJG7IEHkMEybuHAjpEZH335ejDtK79qnYF154QY1Q/Oqrr8RbaNMu0mpDRAgJEmGHM8UFCxaoIl+jyAP210mQdsX6SdiBJk2ayFNPPaV/sQwwFJk7glYn5k3HQtSjS1TrFMV8aIgvX6JNgbh06ZJ8/fXXyYQdrHnsJ5ugBECbeeyoxs4YsdNEq6PnXWTtmt0ZG5u08/Ztkf79pXifPpJkCOQ4YgfhOHXqVLWvY8eOsno1qvU855133pHhw4fL448/LhcMhsqEEJMJu759+8qPP/4oTz/9tDKGrVmzpvowefvtt1ONrpgNM02euG6I2IVnyeLXtYwfP17/8pw+fbqsXLnS6X3ZOGHuOrtZs2bp2+kx/7lTp076dteuXZUBeUJCgp7ut4/WaSe0WtQOETv7TEaKqVgruXLlUvWhx0SkedasIoMG6bcV2rJFjSOr6CRiZ4wS3r59W5544gmnM2fdYdu2bbr1jDPjZEKICYQdQv6fffaZSo9cvHhRhf7btGmjalPqwHAzhDDT5InbV67o22H4YvEj+PJDtEADjRRoqHAEhZ256+zSW9i1a9dOnahq5Qn4bIOogWByJuyAJuzQvYtIs7vCzhgtPPzff3J96FCRRYug+NQ+xKJRcdfRQcQONXhiJ/patWqlBJknGNftDaFICAmSrlgUuzdv3lwGDRoks2fP9sX6SHpgsLEJ82Mq1ijmNCNa1FqNGTPG4f0o7MwdsTMKivQQdoi+TZo0SZ9hjL+9tm3b6rc7E3YpNVC4KuySjRZ79FGRTZtknzU1CwOimdZ5s5cMkzGMpsha88emTZukR48e4gkUdoQEL+yKJTbCzp81dhqYGQpvO2126MiRI5N1KgJtHyZY3HXXXem+TuJbYaeRO3duZU6cXk083377rWS1Rq7ho6eRWsTOkbAzNk84q7EDKOvQ0AyZpWhR+V/hwvKJ4X4viMjc48fFcuhQMmE3btw4vd5vxowZqlYwLSDah4yMBoUdIcEFu2KJhBvTNgEg7AA6IHv37q3bmrz00ks2f3NIj6ErVovWsXnHHOTLly9ZgwKiden5+0X0DMLIHleEnX0DhRb5guDSxKIjjObLRouXg8eOyUvoEs+dW25aI3JV8D6oXFmla43CDs1H//vf/9Q23itaZ7G72PvmUdgRElywK5ZI2LVrASfsAOwWCsPrS0R+//13m05FREa02id2xJoHfI7YR+3SIw1rD5oQ+vTpY7PPlVSs/dQUm6kTKVC1alU9lbrO6mEXFxend6RuKl9ehjZsKFrcOgwRtaZN5ZE1a0T75C1YsKCUL19eP+auXbtsjI579eolPXv2VCdKKUFhR0hw4/WuWERZcLaLOo/UPkCI/8GZfYR1ZFOgCTukWD/++GP9Ov62tC861teZl0AQduC9997T5xgjili6dGmH94MViybK1q5da2MjpKVDUxN2+FuvUKGC2t65c6eaWGEUVDjBuVG6tFQVkR+0nRaLPHf0qPyMk5vs2VUa2ZmwmzdvnkyYMEE++ugj1XmeEsb0MaCwI8TEwg5nsKl1xT788MOqUwsfipUrV1bjolCADPEHUUgCC4hvm1kTASTsQIsWLVSXn9bxh0YdQGEXGp2xQGukSW+ioqJk0aJF8v7778vPP/+cLEWsgRSrJsp27NihizlXPOwcpWNRy7xx40abej2ke/PkySM4Mt4Nu9u3F4tVTDYRkZWok920yUbYoWvfUXoXJ0uI4LkascN1LTpOCDGZsLPHUVcsvoiHDBki33zzjTpjxAdU//791Rnr0qVLvbVu4oupE8BPI8VSApEGRDQAmioQFaGwC42IHUQejIH9BcQcpqCkZudUu3ZtPQKuiShXO2Kd1dkZhR0idlqDECpN19WrJxe++UY0R7sCEGq1a0u+337Tf17GiJ3ekGFtuJg/f77Lwg5C0z6KRwgxobDD2SDmexq7pxyRIUMGVT/y3HPPqa4tEsBTJwIwYgcKFSqkOmM1XnzxRZvCcAo7c4HGGa1m12ixFMhowg6sWbPGK8LOmALVInbGeufDxYvL/Uj/ajtv3pSwTp1kSsaMEmkVcIgeQpht3brV5rlSSsc6EnFMxxISAsJu1KhRqhbE0dk0jDq1Id4ksAkGYQfgy4XUPsDfHXwUtRMHV9JcJHiAiEFktnPnzvLGG29IqAg7mBRjRJmziJ1R2KEsAcINEuwhNFdURfVdEk+eOCF/iAjujRMg+OIhMm8EDRobNmxwKWIHKOwICQFhhw8edFg5Al+2GMkD8WdWzDJSLFiEXWRkpPqy14rUNeBvZr+PBD/4/MBIOXR6BosYRWQZoAYZNWnu1tjh71hrFIFwM47Tw/GNXo2I2GlWJ5jLsv2ll0SmTBGJjlb76kHsicipn3+WLVu2ODRC/vDDDx2ug8KOkOAmzd+I+FBJKQWG6QELFy4Us2KWkWLJauwCVNgBpPTtHfWZhiWBFrXDewqzVo0pTVcidgBOAxqwldKaM2JjYx1G7DSUAO7aVWTFCrmZM6faB6Ogx0aNknDDaDY0sWE2LUAdtP1IMsBULCEhKuwwtNrRmZ2xk83RtAASWARLxE4DtXbG6AeFHQnUdKy7qVj7Ojt7A2SUvSBybR+xA1q0UB54QC4uWSKrrfujEhLkqUWLZAKi3iLKvgU1qgBRxddffz1Zl7xmNq95SAIKO0JCQNihqHnmzJnODxwermrtSGATbMIOXYpGbzvY6xASyMIuJibG5c5eR9YumsDCZyrGq6Uo7ETkrooVpXVsrEwyHOMVWKJERUn+8HB55ZVXVAQQfP7553q9Kjh58qS+jZpWrcyBwo6QEBB2r776qkydOlWmoK7DASjOTa/5jiQ0UrEa8LVbsWKF8kV87LHH/L0cQhTwstNmtRqFHSLMrk7mgQC0N2g2jizT6uyQitXEFlK1Ro89PFepe++Vl0WkC6Jw1v014+MlrGpVyXfkiPLm00AED4bI9mlYPK8WaaSwIyQEhF2VKlVk0qRJaobnI488okaNoYsLYfwffvhBXnvtNX1uIQlcgi1iZ4wYN2vWjKPsSMCANKnW/IDaNW0aj6tpWGfpWGNKVKuzu3Xrll6DZ4zWaaCxC0zHewXr0W5ATd2DD0qX8HA92o3P7cGDB6tt+/Sx9twQkpwkREhwEO5p59ry5cuVV1Lr1q2VYTE+eBBRue+++5LNWiSBR7AKO0ICPR2r4a6wMzZQ2EfsjA0UCQkJToWdcQIFTE2qiMhZbaIH/O66dJHvCxSQrBkyqF0ob4CZvDNhB4ypX0JI4JJUiesBcGSHHxJ86+BuDqGAMWKOioBJ4EFhR0hgCbuUInZGyxON1IQdOCUil+fPl9wffSQyKan6LvbLL2VXiRJy/8GDgtgipgdpqWQthWzfQHH33Xe79VoIIUEo7Iyjf+xnPJLgqrFLiIqSiIgIP6+IkOAFogzlARgtpuGugTbSqBihFxcXl2LETiOlVKwGmiWKly4NA074Bol0764id4UPHhQYNrUQkV9++cVGmAZaxO7ChQtqNnmmTJlk3rx5yi+VEJIcOruGOMaIXWLGjH5eDSHBDQQUylA8idjh5EobpYbGCJuOVxcjdnhOYyeucUybdO4ssny5SL586mpx61iyMvv2yerVq22OYTy2txooUIM9bNiwVMdR2gMxh3njP/30k2qcIoQ4hsIuxDEKO0umTH5eDSHmS8e6K+y0Wa7PP/+8So/CLsXdiB1EnDEdq43j00G6F+bqVVB9J5IVggvjyPbvV9ejoqKUkbEzLzs0hsADz10TehwDUbcRI0bImDFj3HosGjg0HBkrE0KSoLALcYypWAo7QgJD2GFuLOykmjdvbrPfkbBzNnbNmI69//77k98BgnDVKrnYpIm+620RmYMoXt68Shw6E3YYF/nOO+/Ik08+qTz1XAU2WDBGBmjWcIfLly/r22fPnnXrsYSEEhR2oT4r1iDswrJk8fNqCAl+vCHsnOFqKhYgMqalhxs3buz4gBkzSuxPP8m72bJJonXX0yKy4MIFFNVJ3rx59WkXsEXRQMOcZruy3xrlcwXj3Fp3Hgc0rz3gjpgkJNTwSvMEJkxs375dhcoTE7WPhyRatEBZrjlnxeKCs0jNxT0YiY+LE61dgsKOEM8pWrSoiqIhXYg0qjab1RvYR+zQQIDxjo5o1KiR/PPPP6rWztl9QFh4uBx6+ml5fMoUmW1Ny5a9elU1WUTMn69eD46DCz7fMY1i3759DqdVpMbWrVv17X///VcJw+joaJcey4gdIekk7DCOpkOHDg7faAjla15LJDCxWDvvQHhWfKQTQjwBn3tDhgyR/v37S8+ePfWxXN5AmxerpTMRrUvJpNvV6T+Y4AJhBwc9VM2pR506JVK/vvQoV076WMs2kI7NkSOHnMJtVlKaGZ5SxA4i8fDhw1KqVCm3hR0jdoQ4x+NPHMwdRJ0F3tx4oxovFHWBTyKFHSFep1u3bsq4/e23UbXmPSDijFE7Z2lYd2nQoIFqmNglIiguOVwcvbIq1yq9t26VMdYviz179iRLoboq7BDZMwpCcODAgTSlYhmxI8SHwg5v1L59+6paDBJ8hCHlom3TnJgQr+HNSJ2zOjtvCTv45tWrV09tnxeRZQMHirz0kn77qyICg5EDmzcnE3aupmKN0bq01NkxYkeIa3j8yYMCXYwVI8FJ2LVrd65Q2BES8PgiYgeeeOIJfftedNHCzPiTTyTRalr+mIg8/v77csraOOFuxM5YX5eWiJ1R2GEmOTNChPioxg4zBpGKXbVqlTLmRDjfCGpMSOASfv36nSsUdoSErLDD7O/r16+rGjq9279bN7lRpIjcaNpU0H5R6MIF6fTJJ8oSZb2bws5RxC6twg6lPjA49mZjCiFmwWNhN2fOHPntt99UdxYid8ZCXmxT2AU2ETdu3LlCYUdIwOOLVKw28aJPH7RJ2JLpscekcd688tmpU4I2h9gbNwQ5mk4iMjcNqVh8V+C50IzhTirWWGOnpWMp7AjxQSoW7uPDhw9XhcLocDp06JB+OXjwoKeHJz4E8ywjbt68s4PCjpCABzYm2rgxe888X5GpYkV5ALV31uuY0vq1iAxFnfXJk6mmRSHKtOhchQoVlAEzwHdGfHx8qs8PWxTYahlhAwUhPhJ2eMO1bdvWZ4XCxHfgg9JGylHYERLwwJpk9+7dSijlzp07XZ6zbNmyckFEYHP8mWH/CBGZZbHImVTmyG7bts1mCsY999yjtiEI4WfnbrQOsIGCEMd4rMY6duwoc+ciIE+CeZyYgsKOkKAAQsvRFApfoY0nQ2zteRHpjzo3623PoKv28ccxzNWl+rpKlSrpETtX6+yM9XUajNgR4qMaO5xxjR49Wn799VcVYrdvnhg7dqynT0F8xLVr1yjsCCEuCUkj74vI8UyZZIr1MyQLInI1aoj89JNI+fLJHq+NIANVqlSxmTaBOrtHH31UlYbs3LlTzaeFEXNqETsKO0J8JOx27NihD5jGm5IEDxR2hBB3InZGLj/8sNT96Sflb1cQOw4fFqlVS+Sbb0TsZtOuX5/UQ4sRaxUrVlTdt/YRu88//1w6d+4s+fLlU+PLMmXKlGLEjqlYQnwk7JYt08ppSbDBVCwhxBXQfQqbFaOYevDBB2XATz9JdRHZlD+/5IPtCQRY06YikyaJvPCC7jmndb8iCIBonaNU7DcQhNYuW3je1YJItMJULCHpKOxGjED5rGNgdzJ0KPqmSCDCiB0hxJ2o3YoVK9Q25tXWrInJsiLHRWRMixbyAcaFLViA+hyRF18UQVPEW2/ZpGFrIF0roqJymTNnVieXu3btUmnYTZs26fc7duyYzXMzYkdIOgq7+fPn21xH6zqsTvDGv/vuu00r7CZOnKguwex+TmFHCHGnzk4TdiVKlFC1cBr/Inr23Xci/fujsDpp5zvvKHG3qUQJ/X7Vq1fXT/pRa7dy5UrVFQtRd9rQfGEv7FhjR0g6CjtHbuI4u+rUqZO0atVKzMrLL7+sLnitsbGxEoxQ2BFC0lJnh1SqcT64mj4By6sPPhApWlSkd28YZYrMni3Nc+ZUzRYXDRE7AA8+CDttgpERRuwISTs+MZ/Lli2bMi02a7TOLLDGjhDiKhgZaYzeYYIExo+J/fQJTBv6/nuRjBnV1Qrnz8tq2JzkyKEifRpGc2VMMDJy/DgSvHdgjR0hruMzV2FMosCFBC6M2BFCXAXNEh06dFC1dchWgPz58+sRO9TJ6cDXbtkySciJCbMiMEBZcu2ahBkyPFqNnmZ072rETrNKwYmpsbuWEOKlVOyECRNsruPNjTf5F198IU2aNPH08CSdhF1iRISEG7ylCCHECKYLwZLECIQdpmBAYCUrS6lRQ34ZOlRK9+mjZszmxPjCBx9MskN57DHJmTOnivzt2bMn2XOlVGNXrFgx2bdvnx61M9b6EUK8IOzGjRuX7M2PtnhMpBg0aJCnhyfplIpNyJDBd+FbQogpQXerBtKx9vXGSw4flk4i8oOI1MGOq1dFWrTQ7VCQjnUm7BAkQJOFfcQO6VwKO0J8KOzQAUuCP2IHYWc7M4QQQlJGS8UCZGpKly5tc/vGjRvlvIg0RHlOy5YS88MPd+xQDh+W2jVryrRp05IdFxHAixcv6jV8RmEHtwUNYyctISQJj4I0sDZp0KCBbj5JglfYWQwu74QQ4q6ww0n+rFmzbHzrtO+GnPnzSwwaKl599c6DR42S1gsWiLEABN52jtKxRmFnHG/GwAIhXhZ2mAu7fft2Tw5BAqV5gsKOEOJBKrZHjx6qBKd+/fpy6tQpVRenRdTuueeeJDuUMWNEPvooaRsOCj/+KH9ERUk26zEeeeQRh8JOq7HDmLHyhlm0f//9t89fIyHBhsdlVc8++6x89tln3lkNSVduXLly52w5Sxb/LoYQEtQRO5woamlURO0w79VR+lR69ICzvW6HUjc+XpaLSLXCheXhhx9OMWIHK60yZcro+ynsCPFBjd3t27dl+vTp8scffygncWMoHYzVXMhJwJFgSG+EUdgRQjwQdkbQKWu0MFEROyNooFi+XHXHyrlzcr+IrBGRFZGRKQq7rFmzKmNkNGnATmvv3r3ef1GEhLqw27lzp1SuXFlta51KGlpHEwlMLHFx+nZ41qx+XQshJPgoUqSIZM+eXTU6aP8DdLoaP/9tInYaGC+2Zo1I48Zq9FjU0aNS7/XXpYqIbDIIO3THGiN2OC6iduvXr1fjyBApRIqWEOIlYbds2TJPD0ECQNhFZNOqXAghxDUgqH788Uc1Gqxdu3ZKwEGIIWKnGQk7jNhpoIt27VqRRx8V2bFDoi5cUGnZVgZhh9RuYmKiLuyAJuzwXGjQqFixYnq8XEJCo8buyJEjto7jdreRAAaeUlYiGLEjhKSBOnXqyODBg6V48eLKPFiL2B04cCDliJ1GgQIimBkL82KU+4rIzxhhtnNnso5Yo7DTYJ0dIV4WdngzOxrGfO7cOXUbCVzCrMXOCo4TI4R4SLly5dT/cXFx8ueff6pteNFpfnROyZ5d5NdfRVohVieqqeudw4dFxo+3EXaosQMUdoT4UNgZ3cGN4I2NIdEkcAk3zlmksCOEeEnYAW2Oq9M0rD34vpg3T+bnyXNnX58+kvXtt/WrjNgR4sMau759+6r/IeqGDh1qU7yakJCg6h8qVaqU1sOTdIDCjhDiTYzmwS6lYe2JiJAvateWbQsWyJvWXflnzZIZItLVIOxwzMjISOXKQGFHiJeE3ZYtW/SI3Y4dO2wKZbGNYtZXjS7jJOCIxFBuDQo7QogXI3YaLkfsrBQsVEiGY/asiHwSFiZhFouaN4s43q6YGN0cX5sZC8sTNFdgTjkhxANhp3XDdu7cWSZMmKDXPoQKEydOVBdEJ4MRfBBGxcff2UFhRwjxd8RORO6/H652Ip/iCypvXhl/5oxEJiRIU9z2xRcqPSuxsSodC2GHlO/Ro0elaNGiXnsdhAQzHp/ilCxZUubNm5dsP0yL33vvPTErL7/8smrp/+uvvyQYwYehjZSjsCOEeAhSpQULFvQoYte+fXupDo87nECfPCkNEhLkkvW2AgcPijz0kMjp01IaVilWOK+cEC8KuylTptgUsmpgnt/kyZM9PTxJjzmxgMKOEOKDdKy7ETukWefMmaNngVaKyEMicj4iIukOKAOqW1fujY3VH3Po0CFvLJ0QU+CxsDt58qTDsTJ58uSREydOeHp44iMo7Aghvk7HoqkuX758bh8D9XNTp06ViIgI1aDX4NVXJfOmTSKFCiXdYd8+eXLCBClpvf/hw4dVWhZ13ZhT6y6w58Kc2pYtW9qMQiMkJCdPFC5cWNasWZPMsw77CsB4kgQkFHaEEF9H7BCtS+toybZt26p6O9hmYXSZYvVqkYYNRQ4ckIynT8sqEWlsjdi98MILsmLFCvn888+VOX7GjBldfq5+/frpdeOTJk2S3r17p2nNhJgiYte1a1f1JpgxY4aa24cL6uv69OmjbiOBydWrVynsCCE+F3aeUKpUqTuiDqBBYtUqkQoV1NW8ImoEWZZt25SoA2fPnpW5c+e69TwYi6axceNGj9ZMSNBH7Pr376/C2C+99JIewsYZ1muvvSaDBg3yxhqJD2DEjhDiC6pVq6ZGiyE92rp1a+8/AVK7y5eLNG0qsm6dZBeR8bt3y78i8pv1Loi6deoEkxTXMKZfjdZdhIRkxA5hdnS/YqwYRshs27ZNzp8/L2+88YZ3Vkh8AoUdIcQX4MR+586dqlMVHa4+ASPKfv9dNlgbKGCPj5jbE9ab4VbgjmNBvMH6icKOBDtec3TMkiWLOlO79957JcZqIkkCF6ZiCSG+InPmzG7bnKThSeTjxo3lO+tVyDEkYDtbr3/6KZzwXIMRO2ImvCLsVq1aJc8++6zUqlVLjh07pvZ98cUXshqFriTgI3YWFDdzri8hJMgofM890lZEjRwDMESZLiLdRWTr1q0uHwcTlAgxCx4Lu++++04aN26sOpA2b94sN61jqi5duiTvvPOON9ZIfCzsbuMMNY2da4QQ4i/gxoDZP11QZ2fYPwlp2X9RdZe2bAYhIS3s3nrrLWVEDM8hGEtq1K5dWwk9EgTCjtE6QkgQgiYNgHhbHwg6g2nxoLNnRdIQXIiLi/PqGgkJOmGHAcwPPvhgsv2xsbFy8eJFTw9P0qHGLtENvydCCAkU7P1Tf6ldWz6+6647O15/XQSNfE5SrcePH5cbN27Y7KOwIxLqwg6u4gcOHEi2H/V1cA8nQVBjR2FHCAlCYJBvNEAuV768fF2ypPQ33mnkSJHXXksm7pBlwlxb+5GYFHYk2PGKQXGvXr1k/fr16g2GM6DZs2er0S7du6OElQQi1+PiRJNzFnbEEkKCEHSwFtLGjFlnlGfLlk3eF5FXjHccM0akVy+RxER9FyZVAJjqG6GwIxLqBsUDBw6UxMREadCggYoCIS0LuxMIu1desXlrkQAi/tKlO1co7AghQVxnd/ToUX3qBYQd+FhEhr39tuQeMiQpWvfRRyJIu06eLDdSmAdLYUeCHa8YFL/++uvKlBimlDAphlnxSIS/ScCSeOWKvh2eJYtf10IIIWnloYce0suC4KOqCTtwDNMpZs4UCbd+1U2dKtK5s2xav97p8SjsiIR6xM4YEi9btqzaTuvQZ5J+WAwfXuFZs/p1LYQQklYGDx4slStXlgoVKijbLTTuaVy+fFmkQwcRmOY/84xIQoLIrFmSa9s29eV328HxKOxIsOMVg+LPPvtMnSlhlAwu2J42bZo3Dk3SQdhFGM5wCSEkmEDpT8uWLfUOWWPEDsLu0KFDcqxOHZF580Ssllxltm2TedZpFY6EXWqGxXB8gCF/nTp15Ioh+0GIKYQdZsKieaJ58+Yyb948dcF2nz59OC82WISd4QyXEEKCGaOw+/XXX9VoM3TPdl64UM599plYrCMvHxeRbx2IO9SM21ugODLmX7dunaxZs0Z++OEHn7wOQvyWiv3kk09U23i7du30fS1atFBhcTRPjBgxwtOnID4g7No1fTuCqVhCiAmFHQQYhBqYOXOmLF++XNZMnSrZO3SQTCLSHPfBlArMi7WL2iGt64yzMD+2cu7cOR+9EkL8FLGLj4+XqlWrJttfpUoVuX3bUQUDCQTCr1+/c4VdsYQQEwo72G8ZOXz4sHx64IA8BpN2675mIvI9Urpu1Nmp2j0rTMUS0wm79u3bq6idPVOmTJFnUKxKApIIY6qBwo4QYhKMzROO+Pzzz2WFiI24a2on7lITa0ZhZ9wmxDRdsWie+O233+SBBx5Q12FWfOTIEenQoYP07dtXv9/YsWO98XTEC1DYEULMHrFzhGZIvFJEPnj4Yem3dKmawgOhN19EWjFiR0Jd2MG7Dq3m4J9//lH/586dW11wm4bZLFAmTpyoLglonw9CoowGnRR2hJAQEXYaefLkkRyPPy5Nli6VX0QEbp5NrOLu2vnzKT6Wwo6YWtgtW7ZMQpGXX35ZXfAGTy30H2igmDjaWP9IYUcIMbGww0zY69evKyN9jfvvv19NrVhlFXSLDOLuFKZVNGgg4qSBgsKOmLrGDm8WjBIzhrnHjx+vUrMkMMHvCx1hOhR2hBATC7u77rpLKlasaLMPmaYiRYqo7dUi8ihEmvW2vNu2ibRsiS84h89BYUdMLexgDDlr1izdtLF69erywQcfqP2OmipIYAg7GylHYUcIMQmZM2eWcG2EmCHtai/sELErWrSofn2NnbiT3393Ku7YPEFMLew2b94sdevWVdvffvutmteHqB3E3oQJE7yxRuJlKOwIIWYF9dz2UTtnwi579uwqAPHII4+o0WRrreLuVnT0HXHXogU+NG0ey4gdMbWwg0jIajW4Rfq1devW6mwJHbJa9xEJLK5evUphRwgxLY6EXaVKlfTr+M66++671Xa3bt3Ud1e9evXUdYi7z//3P9wp6c5//CHSvLmNuKOwI6YWdhjXsmDBAjl69Kga39KoUSO1//Tp0y53J5H0hRE7QkioCbuyZcvqQYiaNWsmS9dmyYLWiST25Mghc597TnT5tnSpLu5gvG+sK0+LsIMdmDYRg5CAnBX76quvqu6iGjVqqDcMwBkQQt0k8KCwI4SEkrBD80RMTIzMmTNHnnvuOdXgZ49R2MHHbtDChYIwxSWjuGvWTK6cOmXzONzXHZH23nvvqdq+Jk3Qf0tIANqdtGnTRurUqSMnTpywqWFo0KCBtGoFq0cSiMIuu3FHJpseWUIICWrsLagQsQNNmzZVF0cYhd2+ffvk0KFDckhEibvVmTJJFKJ0y5ZJhqeflgwicsNO3DnKUCG6t23bNvXdGBmZ9HU7cOBAPfhx6dKloLPLIiEQsQNomEB0zhjaRndsmTJlvHF44sMau/ioKBG7lAQhhJgtFZsaRmG3YgWGjiWxQURmPP00DqquZ1y9OtlsWWfp2BdeeEHNUn8aj3eA/SxbQrwBv9FDPBV7O8b48UQIIcGPp8LOnvkQYIsX407qOpKo34pIdCrCbsaMGer/7777Tnbs2CG/o8vWwLFjx1JdFyHuQmEX4sIugcKOEGIy0iLsMmbMmKyhQuOvv/4SC2ahL1oktzMgESvSTES+wXhGJ152xikXAFktrblQg8KO+AIKuxBPxSY4GZlDCCFmEHZRUVEu1bHB/85Z1O7cuXOq5k7q1JGVAwaI1hPbUkS+Ro3dhQvJHqPNTtdwNFecwo4EnLCLj49XTRL79+/33oqIz7lmEHYWCjtCiMkwCjlE6yDaXMFe2GkND1rUDvxTqJCK1mnzKFqLSOkRI9ApoX8vTps2TT7++ONUn4/CjgScsMOZ0Pbt2723GpIu3Lp8Wf/FW2h1QggxccTOlTSshrHhD2lZWHlpbNy4Uf2PtOsyEWlh6IwtuHatSPv2Stx98cUX0rVrV33UZkpQ2JGAtDt59tln5bPPPpN3333XOysiPifBWA9CYUcIMRlpFXZfffWVzJs3Tz0G05MQvNC+27QgBixKwB8i8riI/KB1yH79NUJ8MnwZZJ9rsCuWBKSwg0/P9OnT5Y8//pAqVaqoAcxGxo4d6+lTEC9z2/rBBMJT6AQjhJBQEnZ58+aVHj166NctFovkzp1bzp49q/zoxK5R4lcReQLiLiJCIlBD9+WX8kHmzPIUHuvC8zFiRwJS2O3cuVMqV66smzoacbWugaQvlrg4fTtcm4dICCEmoWDBgvo2piKlFXyHVahQQZYuXSqnTp1SF/sO2J9FZE6rVvLsggUqFdvm6lWZAg87F8TdyZMnVXDEWMtHiKd4/Ne0zI2wMwlAYcd5voQQk1GuXDk17hKBh5dfftmjY2FqBIQdgBedI2uTDw4ckJg2baT13LkSYbHI88iMiEj3VI6NUWQQi0YhSkhA2J2sWrVK1drVqlVLDy2jgHT16tXeODzxobCLZMSOEGJChg8froyBPRVNiNhpIB3rSNht3bpVnvr6a3naYlGCDnQTkY9cOD7TsSTghB3eOI0bN1bmjps3b5abN2/qBabvvPOON9ZIvEwYZh5aicxuMzWWEEKIE2GHBgpHwk4D0yjao0HNeh3VeuNSOT488ggJKGH31ltvyeTJk2Xq1Kmqg0ijdu3aSuiRwCP8uubAJBLBiB0hhKSY1o2IiNAjdlonq7bPHhgWd0RmxFpj3ltExqRw/AkTJqiMF4IkjkyMCUl3Ybd371558MEHHRpEXrx40dPDEx8QcUNzX6LdCSGEpESGDBmkdOnSurA7evSo2q5bt67N/d577z19e7aIXPzgA13cwQ1vlPW2++67T5577jn9vosXL5bZs2dLmzZt5IUX0HJBiJ+FXb58+eTAgQPJ9qO+rkSJEp4envgACjtCCHGdhx56KNm+xx57TB599FG13atXL+ndG7G5O2Tv3Vs2GYTaQBEZKSKbNm6UZs0wuyI5c+fOVRYrhPhV2MFhG3/U69evV63hCFPj7AOO3d27p9YTRPxB5K1bd65Q2BFCSIogVWoPassXLlwou3fvVn6t0dHRMmfOHJW6/fTTT9X34anmzeVFw2OGYGLTO+9Idie1zZjjfeTIER++EhIKeGx3MnDgQNWyjZmx165dU2nZmJgYJexeeeUV76ySeJVoCjtCCHGZGjVqSOHChfU0rJZShXgrW7asvu/pp59WF41MmTIpTzt80U7Udg4fLqW7oWfWMRCKRYsW9c0LISGBxxE7/GG//vrrcv78eeUZ9Oeff8qZM2dk5EgEnUmggeLcGGOBLoUdIYSk+j3XsmVLm9SsKwb82iSmSSLS07C/wOTJMsjJY3bt2uXxeklo4xUfO4AwNELQ1atXlywcUxWwIKpqI+Uo7AghJFUGDRokuXLlUt2wcINwBUTsNOBp93bu3Pp1mIENcPAYCjvid2HXoUMHNSv2n3/+8XgxxPdQ2BFCiPsUKFBADh8+LCdOnFBm/O4KO/BNgQIiY+6Yn6CPtp+DVCwhfhV2iNS9++67UrJkSVWDgCLTadOmyf79+z09NPEBKM6lsCOEEPdBNipPnjwu399e2KnU7Kuvirz7rr7vfavXnVHYbdiwQV0I8Yuwg4jbt2+fKiodPXq0+sP/4IMPpEyZMlKoUCFPD0+8DCN2hBCSPjgUduC11+Q9w5xuTKcYlS+f2o6Li1PNGrjMnz9fv8/169fVdUQMCUmXGrscOXKo+gP8j1buyMhIt85sSPpAYUcIIekDRm06E3pfFC4sQw23DTx5Ul62e3yPHj3kypUrahtOE61bt5ZHHnmEXnfEt8Ju8ODBqt4Aog7WJzdu3FD/nzx5UrZs2eLp4YkPhd3t8HARwxg4Qggh3sM4ZtMmYgcD4+zZBS0Ybxpu/1hEjO6v8IX97LPP1PakSZP05grjfNkLFy6opkVE+DCjnRCPfexQX4fI3LBhw9TZRKlSpSQUmDhxoroE22w/1NgVsG7HR0d7/gdACCHEJYwRO4zdBMMxDUhEj95BvhUtXlwGHjqkrqNhw55jx45JbmuHLbp1//rrL7WN7+Hx48enwyshpo7YISoHHzsUetauXVsKFiwo//vf/2TKlCmq9s6svPzyy6rIVXtDBWXELibGz6shhJDQwT5ip/EG/hl0x9nutUOH5Hnr9s2bN5Md57///tO3MfVJY/Pmzd5fNAk6PA7YVKxYUV169uypD0keN26cEj6YSBFsES2zYxR2CRR2hBDil4hdNkPzhOLtt+EgLzJ6tLo6FcbILgg7fM9qwGOPEI+FHYo4EbVbvny5uqxevVouX74sFSpUkHr16nlnlcQndicJdoW9hBBC0idiF44aZytqigUusEGBuPvgA7Uf48hmbN/uMBXrSNgZj0lCF4+FXc6cOVV7NqJ2EHJdu3aVunXrOh1yTPzLjcuXRSvntVDYEUKIX7pkHQoyiLsxY+RKXJxk/fRTtavzpk0i42CI4jhiZ8yKMWJHvCLsvvzySyXkkoWVSUBy29A1ZbHzWCKEEOI7jBE1o2WJTaQtLExuvfWWvPXppzJE29e3rxo/lpSkFZkxY4bqhO3WrRsjdsT7wq5p06Zy8eJFZUq8Z88etQ8zY7t06aJ3/ZDAFHb0sCOEkPRDpVxdsEWJyZBBdcneEpERhvFj0SLKIgV0795dff9S2BF7PP4r2Lhxo9x9992qYeL8+fPqgm3sY4dO4JFw+bK+HZYli1/XQgghoUq/fv10oTd79myb2zJkyKD+Hykik4sW1fePNAg9zQqFwo7Y4/FfQZ8+faRFixbqD+z7779Xl0OHDkmzZs2kd2/jBDwSCCRaXcxBeNasfl0LIYSEasQO89XXrl0rP/30k/oONYLJTZpI+yxPHrkyHG53SQy1Ru8ADIlTqrGLj49XQwTeeOMNOlSEEJHeiNhNnTpV/SHqB42MlAEDBkjVqlU9PTzxMpa4OH07gnWRhBCSbhi/J8EDDzzg9L4xMTFqPizsTs537CiDhg1TkynAAGta9sL58zYRO/tRY5MnT5ZRo0ap7WLFislzzz3n1ddDTBqxQ9PEkSNHku0/evSoZGVEKPC4elXfpLAjhBDfoqVZs2TJIs8884zLj9PSsRjTictEEXnBcDvyYcXff19OHj+u74MQNIKgi8asWbM8eBUkpIRd27ZtVaPE3LlzlZjD5euvv5bnn39e2rVr551VEq8Rdu2avh1FSxpCCPEp+B5cs2aNmlTkjg0YInYAETtNsEGm9c2eXbQYXZ0dO2Ti7dvKyNiRsHOlcYOYD49Tse+//776g+nQoYPcvn1b7/BBxw7myJLAIszwxqewI4QQ34Lvx1q1arn9OC1iB2GHiJ3GqnvukfYbNwrib6iow+gx9NQiyardDylZBFh27Nhhsw4SGngs7KKjo+XDDz9Uefx//vlH7UNHrHF0CgkcIgzCjs0ThBASmGgRO4g1YyQuf/788hUaI0TU//gS72gVd+9aMzIrVqxQM9udce7cOTUhCpm2Vq1aSenSpdPhFZGgEXYaEHL33nuv2uaZQeASYTjzo48dIYQEtrBD5+vSpUtthB2YZxV3c62NFJBx2f/9F62wynLMGdOmTVMTojS++uor2e5gbBkJXrxievPZZ58pUYfQMS7Yxh8PCTwib8Hy0gqFHSGEBCRaKha89dZbyYQdWCAirZGutV5/DBG7tm0l1sG4SARcrl27ZiPqANK1sEUh5sFjYQd/nF69eknz5s1l3rx56oJt+NvhNhJYRFPYEUJI0ETsHM1nRwmUxs8i0hIpW23H/PnSb906FcUzgpQuRpE54tSpU95aNjFDKvaTTz5RLdXGDliYLVaoUEFeeeUVGTHC6JNN/E2U8cyMwo4QQgI+YmckY8aMykoMdXIav4pIMxH5EbeLSMUjR+R7EXnCEM1DSvfXX3HP5Jw8eVIKFSrkk9dBgjBihxCuIyPiKlWq6F2yJDDA7yqjwcySwo4QQoIrYgfB5+i7dQlmt6Mj1tq4iO35uL/1djRLOBvzeeLECS+unAS9sGvfvr2K2tkzZcoUt8wYie9BZ5WNlKOwI4SQoBJ2iNgZ7U+MLEOE7vvv5WYUemRFmojID9YoHjxmjx075hNhh2jgLWOZDzFP8wRMiXG57777VHoWs+769u2rX4h/QeEshR0hhARvKhb74W3njGvVqsnHTZuKNhW8kTVFazQgQ7lUwYIFvSLsvvvuO5XGLVKkiFy8eDHNxyEBVGO3c+dOqVy5strWfOxy586tLrhNgxYo/ufq1asUdoQQEuQRu9QyMzuzZ1eCDhV1GBzZQER+saZnr1onRg0dOlSqVaum19ilhfPnz0ubNm3UdlxcnCxbtkz54pEgF3b4RZLgi9glhoVJuJMPDkIIIYEbsUtN2OHyp4g0FJHf4G8nIvVEZLE1PVuzZk2bztq0Ruy2bduW7LmJSVKxJPiE3a3ISIRR/bwiQggh7kbsMOnJGRBXWg3eX9Zo3XnrbXVEZHXmzFI8Z07JmzevnklLq7Czr/WjsAsMKOxCNBV7y3C2RgghJHi6Ynv27CmjR4+WsmXLJrvdfgTZZqu408xRKl69KtKokUTGxSlPPHD69Ok0rdG+1o/CLjCgsAvRiN1tpmEJISQofewwwrN///6yadMm6datm8NUrJGtIvKwiJzRdmzYINKwoRTNlk3vanUG/PKaNGkiTz75pPoOMcKIXWBCYReiwi6Bwo4QQoIyYmcUebAbQyOEo1SsEUyDndmhg0iePEk7Nm2SL06dklgR1c2aaPQ4NYBpUosXL5Zvv/1W+vXrZ3Ob/fM4s2EhQSTsYHjboEED2b9/v/dWRHzG9bg43awyIZXOKkIIIYEZsUtpn33E7pdfflF2Y8OGDZPemOG+fLnIXXep28pdu6aaK7JZLLJv3z5p3bq1PP7446psR2P9+vX69uTJk23MkV1JxR4+fFi++uorm2OSAO6KjYqKku3bcR5AgoGb57USWhELhR0hhJjC7sS4z1hjlytXLpVGxUWnXDnYWYjUry9y5oxUt9qi1ChbVi5b74JRoO+9955+DCO7du2SihUr6s+VkrBLSEhQlipnz55VKWNHwwxIAKZin332WWVQTAKf24Y6Cm3sDCGEkOARdkabktQidk4976zi7rI1KljD4HkHkHrVuHJFszpOAhMsnEXs7IXeqVOnlKjTon0kSHzsEJadPn26/PHHH2o+bGY709uxY8d6+hTEB8KO5sSEEBJ8qVhHZv9Zs2a1aXbQBFaKnnfly8u0du3k2RkzBInZB6w+d43hmmAYD4YZs86EXWoRO2d1eySIJk8gR2+E0yYCiwTjGzRLFn8uhRBCSBoido4oVaqUTao01YidlZslSyorlKUigpaKmtbIXXdDJC4lYZdajR2FnX/g5IkQItEQUo8wnOERQggJLBxF29Cs6Ihy5cqpQIrFYpGZM2fq+1MTdjly5BAM/tTEXW6ruJt27BhysAgFJkvF/vfffy53xRobLUj6QbuTEMISF6dvh1PYEUJI0ETsZs+eLV9//bXD+8LX7p577nF7/Fj27Bg2JrLD6nOXVA0nUhWp2EcfVeLOk4idMaVLgihip7F79245cuRIsl9kixYtvPUUxFMM7eaRsXAvIoQQEgzC7n//+1+K97/33nuTWY+lFrHThJ0m7rTIneqDXbtWpEkTSbh4Mc01dhR2QSrsDh48KK1atZIdO3booWBjfR3anUngCbsowxuaEEJIYJFatM0ezH61x5VUrJHtBnGnho2tWSOTRARmKdq3xzGkaV3siqWwC9JUbK9evaR48eJq1hzCwSjcXLlypVStWlWWwwiRBAzhhrMpCjtCCDFH8wSAsbA97kTsNLaJyKORkWKx3lYXJscwUjCINy2Aw4idSYXdunXrlJlh7ty5JTw8XF3q1Kkjo0aNUoOKSeAQYXgTUtgRQkjggu9Sd2jUqJF89NFHaaqxs2djQoLsGDdONEv7B0XkZ4O4w9QpfPdjzJgRCjuTCDukWjUPHYi748ePq+2iRYvK3r17PV8h8RrhBmHH5glCCAlcIJ7cAeVPPXr0cCtilzOnSrgmAxG5fVmySEMRuWDdV09EfkKjhogqv6pVq1ayx1HYmUTYoWBz2zYEb0Vq1Kgho0ePljVr1qgoXokSJbyxRuIlooz1EDQoJoSQgKV06dL69nPPPefy41ASldLsVvuxoEuXoqIuOQjMbBGxEXf1reJu+S9IziaHNXYmEXZDhgzRTQgh5g4dOiR169ZVg4cnTJjgjTUSLxFlfJNR2BFCSMACgfbXX3+p71F3JjgZZ7tiCkVqPPTQQ/LwwzA7sUUL2GwWkUdEROuNfUhEfkQ00MGxGLEzSVds48YYQJIEfHT+/vtvOX/+vOq24eSJwCLaGNqnsCOEkIAGTYi4uAOEnWZJ4oqws4/y2Qs7sMkq7n5HXZ7V8w7irjnEnOExWmOF9t1PYWcig2Lk7SnqAo9oows4hR0hhJiOYsWK6dvR0dEuPcZRLZ79iNCNaNAQEW3iOGxRFjqI3BktUCjsgihi17dvX5fv604Imfi2EDeTtUVdQWFHCCGm44MPPpCff/5ZjfMaM2aMS49JrckCpVWYfIELxN1vIhJrrb/7AYMIEK2z3hedua+88orqyHUm7E6ePCnNmzdXXbk//fST29YuxAfCbssWlFSSYOLatWt6q7qCwo4QQkwHmhYxOACRs7vvvttjYdeuXTtp0qSJzJs3T13fgBIsq7jLZk3RQty1tIq7AQMGqPv179/fqbBDcGjjRsQARcaNGycDBw5MwyslXhV2y5YtS8vDiB+5evWqrbBL5QyNEEJIcFKoUCG37u+oxk6jdu3aydK6663i7leruEMUbwFMkg3izpGwQ6PliRMnZM6cOfo+TeCRAGqeQCesM1BnN3ToUE+fgng5YncjIkIyuGl+SQghxJzA9sQZml+d/X3+xIQKq7jLahV6863iTquysxd2SA8/9dRTNvtYhxeAwm7+fPwqbWu5YHkSGRmpwsAUdoEj7DRL4ltRUeLeFEJCCCFmxVmzI2rf7rvvPqfib50hcpfVKvSgCFpZRZwjYbd27VqbfRR2ASjsHNXbXb58WTp16qTcqUngCLt81u14FzulCCGEhC4FCxZUQZqUOmzXGSJ3WUSkiYh8LyJH9+9PJtowxMAeCjvv45N8XLZs2WT48OGM1gVojd1tCjtCCCGpROwKFCigb9tH7N5++2355JNP1ECCPvPmKXEXZ73tMUyuuPdeeW/kSJvH7Ny50+PRaSR1fFZodenSJXUJBt5//30pX768Go/25Zdfihm5Fhd3R9ilMhiaEEIISUnYIYDTrVs3WblypbRu3VrWWKN1mrh7NDFRvkOkL5XnMPrekQBJxdqPDYPrNLpevvjiC9UiHejs2LFDvvrqK9m0aZNaO8arNGvWTPnrmImbF7WBMCIJFHaEEEJcSMVq2KdijdfDw8PVPNvp06eraN0iOGqJSDMR+VZE2iDl6uQ5KOwCUNjBg8YIfsF58uSRjh07yqBBgyTQ2bNnj9SsWVOZKYKKFSvK4sWL5emnnxYzceuCNsZZxEKrE0IIIR5E7OxNhadNmyZ9+vRRzRYQd79YxR3GjsEB70kn4g41+a7y+++/y4IFC6R3795SsmTJNLyi0MDjVCw6YI2Xf/75R/7880955513JGtWrQ8z7SDMC4dq/IGhDgC/VHsmTpyoxqhAnKE4c8MGWCi6BtKvy5cvl4sXL8qFCxfU9rFjx8RsxBsidpYUPIsIIYSEFmXLlnW4P3/+/C4LO3w/o6QJrBSRpqjttt6GyRTf4BgOngPfu66AjtpGjRrJpEmTgiIb6E/Cg6HoH1E0iDdHzJ07V7lYDxs2TDZv3qzu27hxYzl9+rR+n0qVKikBZ385fvy4lCtXTnr27CkPP/ywqhN44IEHJCIiQsxGgvGsiFMnCCGEWHnmmWekZcuW6rvSvo7OlVSsoyaMFdZU7DXr9ZZOxB1q8V1Jx964oQ0tExVAIj5MxTqbG4tfMCJo99xzj/qDyZkzZ5qOD2WekjrHLNquXbtK586d1fXJkyerOXnI9WtjSrZu3Zric7z44ovqAp5//nlThniNwi4sC5rSCSGEEFHBDC0bhvpyrfGxcOHCLkfsjJpAmxG/3CrufsJ0C6t58RcYU4bMkd0YNARZUMb1+eefO/TMs7dFQTDn1Vdf9Upm0Gx4xccOkbKEhAQpXbq02rdv3z71h1KmTBkVNu3Xr5+sXr1a/eK8CX7RaHow1vKhxq9hw4aybh3cdVwD0b277rpL9u7dq9K4EIfOwJmF8ezCnfoAf5Jw5Yq+Hc43AiGEEAf8+OOPKoLXoEEDmwievdhy5muHaVTInOXOnVuaNm0qGEDa3CruUN3dFrXtIjLc8Bhkz3DRNMWoUaPk8cchA+9gH9XD80CAjh8/3guv2lx4nIpFNA5CCr8UiCxc/vvvP3nkkUfU8GDUqz344IOqqNLbnD17VgnKvHnz2uzH9ZMnT7r1GiA6n332WZkxY4ZuyOgI/MHFxsbqF+MZTSBjidOa0CnsCCGEOAa+dEeOHFHfhSkJOWcRu8yZM0uHDh2UMNRYau2MTbRef9N63RF///23Gm5gn2lzZGT84YcfuviqQguPI3ZjxoxRnSrGXDwEz5tvvqkKHXv16iVvvPGG2g5U3InuITpoTD8jYhcU4s4g7CJjY/26FEIIIcGFqxE7Z8IPXbID4BtrvT4T2T0R2e7k8bNmzVIRQwRv4C9rrJsnPhZ2CIXiB26fZj1z5oyepkTO3hdjQxDqRcr31KlTNvtxPV8+bYCWd8Efq7MzlYDmmlbCSmFHCCHEPVytsUuJD0Skgoh0sFqh/CAi1ZB9c3BfCEd0zP7yyy9qRClJ51QsjAnnz5+vUrC4YLtLly56jhx1a6VKlRJvg198lSpVZMmSJfq+xMREdR3edOQO4QZhF50jh1/XQgghJLhwNRVr721nzwvQBNbtYlaPO0cRpvfee081XaJEiqRzxO7TTz9V9XMw9IXPjDpoZKQyKNbMi9FE4egX7ApxcXFy4MAB/Tq88pB7xy+8SJEiKi2K56patapUr15dFVLCIkXrkiVJhBtaxSnsCCGE+DIVCxDg2b9/vxJpGmiBaIVauqxZJeuVK1IftXIi8rJPVh2aeCzssmTJIlOnTlUi7uDBg3rrMvZr2HvjuMPGjRvVmC8Nrb4NYm7mzJnStm1blfZFHR8aJvBcmBxh31AR6kQahF2UycalEUIICcxULJop7EH/69ahQ6XOkCESduuWvCQi20RkitdWG9p4zaAYQq5ChQrqYhR1nlK/fn01w9X+AlGn0aNHD/n3339VO/T69evV9AliS6ShVZxdsYQQQtzBFYNid7DUqCFXrX534GMRaZiW41iMjngpgzIxBH/sO37NhscRO4CaNlzQRIEaNyMwCib+J9LYvMLJE4QQQtIhYhcfH+9wP6Y8hdWsKR/06CH9rBMpvhWROiKy0411YXzp66+/7tJ9W7durf5HX4CZy7U8jtgNHz5cWZlA2MFXDl0sxgsJDKKNby4KO0IIIT6usXMk7NDYuHPnTvV4HPONmBjRJsDDr+FnzKh1Y11Dhgyxub5o0SIVmXMnkmc2PI7YYUoD0qLt27f3zoqIT4i2NrYoKOwIIYSkQ8TO3ups7dq1NtczZ8smz5w5o8aPwfqkCKZfiEg9zIp3c4049mOPPaa2MSINrh2hiMcRO/zSatWq5Z3VEJ+AM5cMFHaEEELSiH2EDh6y7go7jPy0B7Ner1nHjh227qsiInPSIFDefvttfbt///7qfxgcz5s3T/744w+b+9qXjZkJj4Xd888/L1999ZV3VkN8AkLhGMCsQ2FHCCHEg4idqxiFnaP0rTa1CmMGmorIRet+CL1P3HieixcvKjNjDS0VC1H31FNPqTGnRnwxNME0qdgbN27IlClTlBpGR6z9L3+soeuF+Idr164pl29wKzxcol080yKEEEJ8KeyMLhq70eAgIr9amylgZnwGdXQuPE+/fmjBuIMWkevWrZvD+2PGPXxyf/jhBylatKiYCY+F3fbt23WfOhREGgkLC/P08MQLwLBZE3Y3IyPFsyZ1QgghoUZa7U2Mws5RXV7GjBltri+zjhybbU0pvm4VdzAxTonPP//c5rrFGrFzpkPWrFmj/u/evbtNpM8MeCzsli3Dr4EES8QOwo4udoQQQtIjYpc//50eV0yhSk3Yga9FJJfV2w6MF5FzIvJlCs+DWjpHwu7iRS2565jVq1eL2fCKjx3YvXu3HDlyxEadQyk3b45MOfG3sLvLuh3voakkIYSQ0COtwg52JKhzw/eQI2NgR8IOTLSKu+HW63jkBasdij29e/dOts9iscjXX0MipkymTDYV6KbAY2GHMWKtWrWSHTt2KCFnH/60V9Ek/bkaF6dH7G5T2BFCCEmnVGz27Nnln3/+UTVvjo6RIUMGp48dISK5ReQVq1iZJyKPishKu/t9+GHyRK3FYpF27dqlaeSZhHpXbK9evaR48eJq6gSU765du2TlypVStWpVWb4czjTE31y/dElX8LdTeBMRQgghjnDV3sQRkajtdiIMnUXsNHqJiOa7kdEasXPFYO3KlSsurc2METuPhd26detkxIgRkjt3buVRg0udOnVk1KhR0rNnT++sknjELcMEkAQKO0IIIWmgSZMm6v/Bgwd77ZgpCTtYlOw/cEDWd+smP1n3oYd2kYhUT+W4F1ycfGVGYedxKhapVhgMAoi748ePS+nSpVX78N69e72xRuIh8YbiUQuFHSGEkDSwcOFC2b9/v8MmCF8IO3jc3X333ZIzXz5pIyI/iEhj7LdaojQUkU0ePn8mEwo7jyN29957r2zbtk1t16hRQ0aPHq3aiBHFK1GihDfWSLwYsbOYsJ6AEEKI70FKtWzZsl61Mkupxk4zL46NjZWbIvK4iCyx3pZdRH4TkSSztbSzfPlyOXnyZLL9ly5dkps38awhKOzQ8aIZAULMHTp0SOrWrat8YSZMmOCNNRIPuX3p0p0rJjw7IYQQEpykFLHLly+fLuzADRFpISIrrLfnFJHfReR+D9fQoQOc8+6AYFWBAgWkSJEiqdqlmFLYNW7cWFq3hle0yD333CN///23nD17VjVTPPzww95YI/GQhMuX9e0wg8s3IYQQEkjCThv9hXr9Ll26JDM2vmYdPaa5z6Frdikyhh6s4fffIQ/v8Oyzzyp7FugY9AuEnLBzRM6cOTl1IoBINHQHUdgRQggJVGE3adIkGTZsmJpkhfo6bXSpkasi8phB3CEtC2n2oJfWdPToUX0b4i7Y8ImwI4GFJS5O346w1iwQQggh/sa+xg6ZvzfffFPV8mm0bNlSF4CzZ2PYmMgVayPFH9b7oIVzsYg08rK1SzB68VLYhQAUdoQQQgKR1HzsQK5cudR0K3jkGk2HkZZtJqJboeBIC611eO5iFHAUdiTgCbuGP/8kIrMjaE0IIYT4n5S6Yo0UK1ZMNWbal3mhb7W1dSoFiBGR70SkvZvrMI5DRX2fhtYcijo8+Ph99913+j7Tz4olwSHsoinsCCGEBPkMWiPxIoI43nV0uEqSsJmFrloRGeOGsNOih44ido0aJSV5Fy9eLOXKlZNNmza5LErTG0bsQoDw6/hzTyI6Rw6/roUQQgjxxqgyI5BfnUTkI8O+0SIyDsENFx6vedYhGodBC/pxExJkyRLNPS8JpIWnTp0qgQqFXQgQYTBZjMkJ5x9CCCEkMEyPvYVFRDDIdJBhX28RmYOghoup2M6dO9vsh7Br2BAzLmw5ceKEBCoUdiFApKFVnBE7QgghZhR2Gu9ao3e3rdfbYsJEhgxqFFlqEbtZs5DEvYOzerpArrOjsAsBogxFofSxI4QQEihgxry7PPXUU6ne53Nrdyw870DNGzdkjYgUdXL/PXv2yODBg5Ptd9YVS2FH/Ep0PEpLrXBWLCGEkAChZMmS0r17dylatKia2+oKMDHGIITUWCQiav5Vrlzq+r0iskFEHnBw3+bNmzucMoFpWo6wWJD4DUwo7EKAmNtaQJrCjhBCSGABoXb48GGpV6+eS/eHrx1q3D788ENp0SJl17oN+GfdOjlk7b69S0SWWbtoXeHgwYMO9zNiR/wKhR0hhBAzER0dLT179pT27V1wrCtZUtoVLy5abytMSr4SkeEudsw6gsKO+A2EizNY/wCVvItOrTeIEEIICX7Cw8OVoTC4Ehkpj4rIp4bb3xCRuYh3mEzY0aDY5KCFW/ujvRERIVnsXLsJIYSQYMVZrVt8fLxcuHBB8uTJo4s8BDe6oW5ORD6wRraeFBFMpW0lIgfceN5AHjXGiJ3JuXr1qi7sbvqgrZwQQggJRBuVPFZRZz8mbDyaJUTkkvU6mir+EpGmbhz/hsFGLNCgsDM5165d04XdLQo7QgghIYj9jNlfRKSaiOyyXsewzZ9EZJiLdXf4bg1UKOxCSdixvo4QQoiJcNV2JNwQsdPYb7U+mWfY96aILBSR1Kz8N27cGLCWJxR2Jufa5csSY92+TWFHCCHERBjTre4KOxAHw2MRGWCdNwuaicgWmBqLc/bv3y+bNm2SQITCzuTcOHdO376dAU3ehBBCiDmoX7++MheGwOvdu7fUrl1bFi5EzC3lVKw9Y0RU1+xZ63VMqFhpFXzOHrloESyQAw8WXZmcWxcu6NsJFHaEEEJMBAQbhBzsR5xF5UBKt2n8ISKVrB53D1oF0nsi8pCIdBCRM3b33717t/r/9u3bqgs3Y8aMEggwYhdCwi4xQP7oCCGEEG+SmnALc9Hq65h1DNlIfGda9yGSt80q8Iz8+++/cuXKFTUWLX/+/LrQ8zcUdiEk7CwUdoQQQkKQMDc8XBOs5sWPiMhJ6778IrLU6n8XYxB2Y8eOVePQLl26JE8//bQEAhR2Jifh8uU7VzhOjBBCSAgSnkpEr1o1mJ/YAiFXUUR+M+zrKyKbrCnb48eP29TZ/f03rI/9D4VdCAm7sCxZ/LoWQgghJNAiduvXr5dGjRo5vO20NRXbB82I1n3lRWSDiAyC7cn69QE3jYLCzuQkXrmib1PYEUIICUXCU4jYVa1aVW7evOnwtkyZMonFOq2iitUGBUSJyDvWztm7A2x+LIVdCAm7iGzZ/LoWQgghJJCE3RNPPKFuczYibMqUKfo2WiNqiMjbBs+7WiKyXUR6B5CgCpR1EB9huXpV36awI4QQEoqEOUnFzpuXNHfCWcTumWeekaZN70yRjReRISJSV0T+se7LJCLjrNE7OXxY/A2FnckJMwi7yNhYv66FEEIICZSIXaNGjXTB50zYgWwOgiLrrI0VEwz7SuGfTJB5/oXCzuSEGQYVR2XHmGNCCCEktOjSpYvN9TJlytikWR2lYseMGZNitA9hk15WM+MDIvIKdt51l/gbCjuTE379ur4dkzOnX9dCCCGE+IOnn35aec6NGjVKTYmAmXDRohgcdiflarQ+iYuLk1dffdXhsYYNG2ZzfZWI3CsicyUw4EgxkxNhCC9H58jh17UQQggh/krF9ukD0xLHYN7s+PHjlaAbMGCAREWh7zU5mEmbwcF4TueJ3PSHws7kRBqEXYZcufy6FkIIISQQCQsLk169kFhNnfPnz0sgw1SsyYlixI4QQgjxGrdv35ZAhsLO5ETFozk7CRoUE0IIIe5h3zzRs2dPiYiIkECFws7kxBjPLDgrlhBCCEkzFotFihUrppovJk2aZHNbixYtJBBgjV2ICDsMOgnPmNHfyyGEEEKCijAHdielSpWSW7du2ewbOXKkBAKM2JmcjNahxDdgzpjCEGRCCCGEJMeYdo2MvBMPy27nDRsdHS2BAIWdyUPGGS0YXyxyPYUByIQQQghxzIgRI/TJFXPmzNH357JzmnBmkZLeMBVrYjAiRauqu2k4yyCEEEKIaxQpUkT27t0rV65ckfvvv1/fn9GuvInCjvicq1ev6sLuFoUdIYQQkibuueeeVO/DVCzxOdfi4kQbR3wrQP7gCCGEEDMSGSABFAo7E3P93Dl9O57CjhBCCPEqsbGxTlOz/oLCzsTcMAi7hJgYv66FEEIIMRsrV66Uli1byowZMyRzgHjFBkbckPiEmxcu6NsJDoYWE0IIISTtVKhQQRYsWCCBBCN2JibeIOwSAyRETAghhBDfQWFnYuIvXtS3LZm0NgpCCCGEmBUKuxARdpwTSwghhJgfCjsTk3D58p0rFHaEEEKI6aGwMzGJV67o2+FZs/p1LYQQQgjxPRR2JobCjhBCCAktKOzMzNWr+maEwUSREEIIIeaEws7MXLumb0ZS2BFCCCGmh8LOxIQbhF0UhR0hhBBieijsTEzY9ev6dkzOnH5dCyGEEEJ8D4WdiYm4cUPfjs6Rw69rIYQQQojvobAzMZE3b+rbGXLl8utaCCGEEOJ7KOxMTBSFHSGEEBJSUNiZmKj4eH2bqVhCCCHE/FDYmZiY27f17TCOFCOEEEJMD4VdCAg71RsbEeHv5RBCCCHEx1DYmZgMCQnq/+vh/DUTQgghoQC/8U1MpsRE9f8NRusIIYSQkIDCzqRYLBbJZN2+SWFHCCGEhAQUdibl+rVrorVL3IyK8vNqCCGEEJIeUNiZlGvnz+u/3HgKO0IIISQkoLAzKTfOndO346Oj/boWQgghhKQPFHYm5frZs/r27ZgYv66FEEIIIekDhZ1JuXXhgr6dkDGjX9dCCCGEkPSBwi4EhF1ihgx+XQshhBBC0gcKO5MSf/Givm3JpBmfEEIIIcTMUNiZlNuXLt25wjmxhBBCSEhAYWdSbl++fOcKhR0hhBASElDYmZREg7ALy5rVr2shhBBCSPpAYWdSEq9c0bcjKOwIIYSQkIDCzqxcvapvRmTL5telEEIIISR9oLALAWEXGRvr16UQQgghJH2gsDMpYdev69tR2bP7dS2EEEIISR8o7ERk4sSJUqxYMcmQIYPUqFFDNmzYIMFO+LVr+nZ0jhx+XQshhBBC0oeQF3Zz586Vvn37yrBhw2Tz5s1SsWJFady4sZw+fVqCmYgbN/TtmJw5/boWQgghhKQPkRLijB07Vrp27SqdO3dW1ydPniw///yzTJ8+XQYOHOjv5cnOadPk4saNbj8uh0GYUtgRQgghoUFIC7tbt27Jpk2bZNCgQfq+8PBwadiwoaxbt87hY27evKkuGpeNRsA+4PzYsfLgnj0eHSNj7txeWw8hhBBCApeQTsWePXtWEhISJG/evDb7cf3kyZMOHzNq1CiJjY3VL4ULF5ZA5lBUlGQtWNDfyyCEEEJIOhDSEbu0gOgeavKMETtfirtsL74oK5xED1MjPCZGSvfrJ+GR/DUTQgghoUBIf+Pnzp1bIiIi5NSpUzb7cT1fvnwOHxMTE6Mu6UWlXr1EcCGEEEIISYWQTsVGR0dLlSpVZMmSJfq+xMREdb1mzZp+XRshhBBCiLuEdMQOIK3asWNHqVq1qlSvXl3Gjx8vV69e1btkCSGEEEKChZAXdm3btpUzZ87IG2+8oRomKlWqJIsXL07WUEEIIYQQEuiEWSwWi78XEcygeQLdsZcuXZJs2bL5ezmEEEIICWGtEdI1doQQQgghZoLCjhBCCCHEJFDYEUIIIYSYBAo7QgghhBCTQGFHCCGEEGISKOwIIYQQQkwChR0hhBBCiEmgsCOEEEIIMQkUdoQQQgghJoHCjhBCCCHEJIT8rFhP0SayYdwHIYQQQoi30TSGK1NgKew85MqVK+r/woUL+3sphBBCCDG55sDM2JQIs7gi/4hTEhMT5fjx45I1a1YJCwvzmVKHcDx69Giqw39J6MK/E/eoVq2a/PXXXxIs+HO96fHc3nwObxzLk2O4+1h37s/3eWi+zy0WixJ1BQoUkPDwlKvoGLHzEPyACxUqlC7PhTcx38gkNfh34hoRERFB9XPy53rT47m9+RzeOJYnx3D3sWl5Lr7PQ+99HptKpE6DzROEkJDk5ZdflmDCn+tNj+f25nN441ieHMPdxwbb32Iw8XKQ/Wy9sV6mYoMAhN6h1C9duhRUZx4kfeHfCSHmh+9zkhqM2AUBMTExMmzYMPU/Ic7g3wkh5ofvc5IajNgRQgghhJgERuwIIYQQQkwChR0hhBBCiEmgsCOEEEIIMQkUdoQQQgghJoHCLsCZOHGiFCtWTDJkyCA1atSQDRs2+HtJJEhp1aqV5MiRQ9q0aePvpRBCvAwmUdSvX1/KlSsnFSpUkHnz5vl7ScRPUNgFMHPnzpW+ffuq1vbNmzdLxYoVpXHjxnL69Gl/L40EIb169ZJZs2b5exmEEB8QGRkp48ePl927d8tvv/0mvXv3lqtXr/p7WcQPUNgFMGPHjpWuXbtK586d1VnY5MmTJVOmTDJ9+nR/L40EITibx0xjQoj5yJ8/v1SqVElt58uXT3Lnzi3nz5/397KIH6CwC1Bu3bolmzZtkoYNG9rMpcX1devW+XVtJP1ZuXKlNG/eXA2ADgsLkwULFiS7D9P2hAQv3nyP47sjISFBChcunA4rJ4EGhV2AcvbsWfXGzJs3r81+XD958qTf1kX8A1IqSMXjg90RTNsTEtx46z2OKF2HDh1kypQp6bRyEmhQ2BESBDRp0kTeeust1QDhCKbtCQluvPEev3nzpjz++OMycOBAqVWrVjqungQSFHYBCuojIiIi5NSpUzb7cR31E4RoMG1PiLlx5T2O6aCdOnWShx9+WNq3b+/H1RJ/Q2EXoERHR0uVKlVkyZIl+r7ExER1vWbNmn5dGwnOtD2+BJ588kn55ZdfpFChQhR9hJjoPb5mzRqVrkVtHpoocNmxY4efVkz8SaRfn52kCOopOnbsKFWrVpXq1aurVnbUYSAUT4i7/PHHH/5eAiHER9SpU0ed/BNCYRfAtG3bVs6cOSNvvPGGOivDGdjixYuTnbWR0IZpe0LMDd/jxB2Yig1wevToIf/++68qil2/fr1qcSfECNP2hJgbvseJOzBiR0gQEBcXJwcOHNCvHzp0SLZu3So5c+aUIkWKMG1PSJDD9zjxFmEWtNIQQgKa5cuXy0MPPZRsPz7oZ86cqbY//vhjGTNmjJ62nzBhAiO8hAQJfI8Tb0FhRwghhBBiElhjRwghhBBiEijsCCGEEEJMAoUdIYQQQohJoLAjhBBCCDEJFHaEEEIIISaBwo4QQgghxCRQ2BFCCCGEmAQKO0IIIYQQk0BhRwghhBBiEijsCCGEEEJMAoUdIYQQQohJoLAjhAQt9evXl969e0uwE+yvI73X/+qrr8rjjz+ebs9HSDBBYUcIISSo2Lp1q1SqVMnfyyAkIKGwI4SENLdu3fL3Eoibv5dt27ZR2BHiBAo7QojbLF68WOrUqSPZs2eXXLlySbNmzeSff/6xSc317NlTBgwYIDlz5pR8+fLJm2++aXOMK1euyDPPPCOZM2eW/Pnzy7hx45Kl9IoVKybjx4+3eRy+0O2P5eq6tLX16NFDPU/u3LmlcePGDo+F+73yyivqfjly5JC8efPK1KlT5erVq9K5c2fJmjWr3HPPPbJo0SL9MTdv3lSv+6677pIMGTKotfz11182x8XjO3ToIFmyZFGv+4MPPkj23ImJiTJq1CgpXry4ZMyYUSpWrCjffvutk9+Gaz9zV36WaXnNGrdv31Y/19jYWPVzHTp0qFgsFpdfk6u/l//++0/Onj2rHk8ISQ6FHSHEbfBF37dvX9m4caMsWbJEwsPDpVWrVurLW+Pzzz9Xom39+vUyevRoGTFihPz+++/67Xj8mjVrZOHChWr/qlWrZPPmzT5fl7a26Oho9fyTJ092ejzcDyJjw4YNSvB0795dnnzySalVq5Zaa6NGjaR9+/Zy7do1dX+Iqu+++049DrdDBEGgnD9/Xj9m//79ZcWKFfLDDz/Ib7/9JsuXL0/2uiGAZs2apda2a9cu6dOnjzz77LPqcSmR2s/cFdx9zcbHRUZGqsd9+OGHMnbsWJk2bZpbr8mV3wvSsBCPEIiEEAdYCCHEQ86cOYPQjGXHjh3qer169Sx16tSxuU+1atUsr732mtq+fPmyJSoqyjJv3jz99osXL1oyZcpk6dWrl76vaNGilnHjxtkcp2LFipZhw4bpz2O8f2rr0h5z//33p/qa7F/D7du3LZkzZ7a0b99e33fixAl1/HXr1lni4uLUa5o9e7Z++61btywFChSwjB49Wl2/cuWKJTo62vLNN9/o9zl37pwlY8aM+uu4ceOG+jmsXbvWZj1dunSxtGvXzuX12v/MU/tZpuU1Gx9XtmxZS2Jior4Pz4t9rr4mV38vI0eOtDz44IOp3o+QUIURO0KI2+zfv1/atWsnJUqUkGzZsqk0Hzhy5Ih+nwoVKtg8BmnH06dPq+2DBw9KfHy8VK9eXb8dUZjSpUv7fF2gSpUqLh3P+BoiIiJUeve+++7T9yFVCfC6kPLFa6pdu7Z+e1RUlHqNe/bsUddxH9SO1ahRQ78P0qbG133gwAEVDXvkkUdUula7INpln1ZOab32P3NXcec1G3nggQckLCxMv16zZk31+0hISHD5Nbnye0HEjmlYQpwTmcJthBDikObNm0vRokVV/VWBAgVUqvPee++1KXiHqDGCL337lGhqIJVqrNMCEE+erAsgXekKjl6DcZ8mZNx9XSkRFxen/v/555+lYMGCNrfFxMS4vV5tba7+LH3xml19Ta78XiDsHnvsMZefm5BQgxE7QohbnDt3Tvbu3StDhgyRBg0aSNmyZeXChQtuHQMRNYgFY2PBpUuXZN++fTb3y5Mnj5w4cUK/fvnyZTl06JDP1uUJd999t14fZhROeI3lypXT74PXjRo4DazR+LpxX4gdRBlRo2e8FC5cOM3rc+dnmRaMrwn8+eefUrJkSRX189ZrQsMNor3siCXEOYzYEULcAt2SSM9NmTJFpfrwZT1w4EC3joHuyo4dO6pGAqQi0UU6bNgwFVUypvMefvhhmTlzporEodP1jTfeUELBV+vyBESb0GigvaYiRYqoBgakILt06aLug/QjtnEfrBWv+/XXX1ev2/izgQEvmgsQFUNnLUQvBCPSy/i5pQV3fpZpAT9vNK68+OKLqsnio48+0jt+vfWaYHOCNZcvX95r6ybEbFDYEULcAiLk66+/VtYaSHOiPmzChAnKrsId0DXZrVs3ZUmCL3d0lB49elTZhGgMGjRIRZVwH9TgjRw50mmUyVvr8oR3331XCRd0jSK6VLVqVfn111+V6NQYM2aMSk1CYEHw9OvXT4kcI3idiLChkxQRKgixypUry+DBg9O8Nnd+lmkBFi7Xr19XNYUQX7169ZIXXnjBq68JadgyZcqkmpImJJQJQweFvxdBCCGwKkH9FaI8WoSLEEKIezBiRwjxC1u2bJG///5bRXgQsYLnGmjZsqW/l0YIIUELhR0hxG+8//77quEBTQewuoBJMcxxCSGEpA2mYgkhhBBCTALtTgghhBBCTAKFHSGEEEKISaCwI4QQQggxCRR2hBBCCCEmgcKOEEIIIcQkUNgRQgghhJgECjtCCCGEEJNAYUcIIYQQYhIo7AghhBBCTAKFHSGEEEKISaCwI4QQQggRc/B/px0spRi3VXgAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# compute the expected number of galaxies in each pixel\n", "nbar = ARCMIN2_SPHERE / npix * n_arcmin2\n", "\n", "# normalise the maps by the expected number of galaxies in each pixel\n", "she /= nbar\n", "num /= nbar\n", "\n", "# get the angular power spectra from the galaxy shears\n", "cls = hp.anafast([num, she.real, she.imag], pol=True, lmax=lmax, use_pixel_weights=True)\n", "\n", "# get the theory cls from CAMB\n", "pars.NonLinear = \"NonLinear_both\"\n", "pars.Want_CMB = False\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", "# factor transforming convergence to shear\n", "ell = np.arange(lmax + 1)\n", "fl = (ell + 2) * (ell + 1) * ell * (ell - 1) / np.clip(ell**2 * (ell + 1) ** 2, 1, None)\n", "\n", "# the noise level from discrete observations with shape noise\n", "ntot = nbar * npix\n", "nl = 4 * np.pi / ntot * sigma_e**2 * (ell >= 2)\n", "\n", "# mixing matrix for uniform distribution of points\n", "mm = (1 - 1 / ntot) * np.eye(lmax + 1, lmax + 1) + (ell + 1 / 2) / ntot\n", "mm[:2, :] = mm[:, :2] = 0\n", "\n", "# the shear pixel window function for HEALPix\n", "_, pw = hp.pixwin(nside, lmax=lmax, pol=True)\n", "\n", "# plot the realised and expected cls\n", "plt.plot(ell, cls[1] - nl, \"-k\", lw=2, label=\"simulation\")\n", "plt.plot(ell, pw**2 * mm @ (fl * theory_cls[\"W1xW1\"]), \"-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(\"angular mode number $l$\")\n", "plt.ylabel(\"angular power spectrum $C_l^{EE}$\")\n", "plt.legend(frameon=False)\n", "plt.tight_layout()\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 }