{ "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", "ws = 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, ws)\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": [ "# compute Gaussian spectra for lognormal fields from discretised spectra\n", "gls = glass.lognormal_gls(cls)\n", "\n", "# generator for lognormal matter fields\n", "matter = glass.generate_lognormal(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, ws)" ] }, { "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, ws[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": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAG5CAYAAAB1OMuOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACJaklEQVR4nO3dd3zN1//A8dfNlsgUESFmrIhN7E3Vrv6qqrVatFVb9YtqaWkptVullFZ1qSodNqVF7cSIHStGECshRMb9/P649356782QXDe5N/J+Ph559LPv+6aS+84573OORlEUBSGEEEIIoXKwdQBCCCGEEPZGEiQhhBBCCDOSIAkhhBBCmJEESQghhBDCjCRIQgghhBBmJEESQgghhDAjCZIQQgghhBlJkIQQQgghzDjZOoD8SqvVcvXqVTw9PdFoNLYORwghhBDZoCgK9+7dIygoCAeHzNuJJEGy0NWrVwkODrZ1GEIIIYSwwKVLlyhZsmSm5yVBspCnpyeg+wZ7eXnZOBohhBBCZEdCQgLBwcHq53hmJEHKofnz5zN//nzS0tIA8PLykgRJCCGEyGceVx6jkcVqLZOQkIC3tzfx8fGSIAkhhBD5RHY/v2UUmxBCCCGEGUmQhBBCCCHMSIIkhBBCCGFGEiQhhBBCCDOSIAkhhBBCmJEESQghhBDCjCRIQgghhBBmJEESQgghhDAjCZIQQghhh/r168dzzz2X66/zwQcfULNmTbt5jr2QBEkIIYSwQ3PnzuWbb76xdRgZ0mg0rFmzxuTY6NGj2bp1q20CygWyFpvItrS0NBwdHW0dhhBCFAje3t62DiFHChcuTOHChW0dhtVIC5J4LK1Wy8iRIylUqBBjxoyxdThCCPFU+eWXX6hWrRqFChWiSJEitGnThsTExHRdbC1atGDo0KGMGDECX19fihUrxuLFi0lMTOTVV1/F09OTkJAQ1q9fr97zzTff4OPjY/J6a9asyXKh1v3799O2bVv8/f3x9vamefPmREREqOfLlCkDQLdu3dBoNOq+eRebVqtl0qRJlCxZEldXV2rWrMmGDRvU8xcuXECj0fDrr7/SsmVL3N3dqVGjBrt37875NzEXSIIksqQoCiNHjmTOnDmkpKTw6aefcvLkSVuHJYQQT4XY2Fh69uzJa6+9xokTJ9i+fTvPP/88ma0jv2zZMvz9/dm3bx9Dhw5l0KBBdO/enUaNGhEREcEzzzxD7969efDggcUx3bt3j759+7Jz50727NlDhQoV6NChA/fu3QN0CRTA119/TWxsrLpvbu7cucycOZMZM2Zw5MgR2rVrR5cuXThz5ozJdePHj2f06NEcOnSIihUr0rNnT1JTUy2O31qkiy2H5s+fz/z580lLS7N1KHnivffeY968eeq+oihMnz6dpUuX2jAqIYR4vLp163Lt2rU8f93AwEAOHDiQrWtjY2NJTU3l+eefp3Tp0gBUq1Yt0+tr1KjBe++9B8C4ceP45JNP8Pf3Z+DAgQBMmDCBBQsWcOTIERo0aGBR/K1atTLZX7RoET4+Pvz999906tSJokWLAuDj40NgYGCmz5kxYwZjxozhpZdeAmDatGls27aNOXPmMH/+fPW60aNH07FjRwA+/PBDqlatSnR0NJUrV7YofmuRBCmHBg8ezODBg0lISMh3/cM5NWXKFKZMmaLuFypUiIcPH7J8+XI+/PBDgoODbRidEEJk7dq1a1y5csXWYWSpRo0atG7dmmrVqtGuXTueeeYZXnjhBXx9fTO8vnr16uq2o6MjRYoUMUmoihUrBsCNGzcsjun69eu89957bN++nRs3bpCWlsaDBw+IiYnJ9jMSEhK4evUqjRs3NjneuHFjDh8+bHLM+D0VL15cjV8SJGGX5s2bx/jx49X9zz//nNjYWD7++GNSU1OZOXMmc+bMsV2AQgjxGFm1btjL6zo6OrJ582b+/fdfNm3axGeffcb48ePZu3dvhtc7Ozub7Gs0GpNjhtoirVYLgIODQ7ruupSUlCxj6tu3L7du3WLu3LmULl0aV1dXGjZsSHJycrbfV05kFb8tSYIk0lmyZAnDhw9X96dNm8bgwYOJi4tj1qxZPHz4kMWLF/Pee+/h7+9vw0iFECJz2e3msjWNRkPjxo1p3LgxEyZMoHTp0qxevdoqzy5atCj37t0jMTERDw8PAA4dOpTlPbt27eKLL76gQ4cOAFy6dImbN2+aXOPs7JxlqYmXlxdBQUHs2rWL5s2bmzw7PDzcwneTt6RIW5j48ccf1b5s0PVn/+9//wN0P2gDBgwA4MGDBya1SUIIIXJu7969TJkyhQMHDhATE8Ovv/5KXFwcVapUscrz69evj7u7O++++y5nz57lhx9+eOzcShUqVGD58uWcOHGCvXv38sorr1CoUCGTa8qUKcPWrVu5du0ad+7cyfA577zzDtOmTWPFihWcOnWKsWPHcujQIZM/wO2ZJEhC9dtvv9G7d2+1OXbUqFF88MEHJteMHj0aJyddw+Nnn32mjmoQQgiRc15eXvzzzz906NCBihUr8t577zFz5kzat29vlef7+fnx3XffsW7dOqpVq8aPP/6Y7ve6uSVLlnDnzh1q165N7969GTZsGAEBASbXzJw5k82bNxMcHEytWrUyfM6wYcMYNWoUb7/9NtWqVWPDhg38/vvvVKhQwSrvLbdplMzGEoosGYq04+Pj8fLysnU4T2zz5s106tRJ7WN+4403WLBgQYZzZbz66qvqXyCffvopo0ePzstQhRBCCItl9/NbWpAEO3bsoGvXrmpy1KtXL7744otMJxIbM2aMem7WrFk8evQoz2IVQggh8oIkSAXcgQMH6NixIw8fPgR0M6N+/fXXODhk/k+jcuXKdOvWDdDN4bFs2bI8iVUIIYTIK5IgFWBHjx6lXbt2ah3Rs88+y48//qjWGGVl3Lhx6va0adPsYtZTIYQQwlokQSqgTp8+Tdu2bbl9+zYAzZs3Z9WqVbi6umbr/rp169KmTRsAzp07x8qVK3MtViGEECKvSYJUAF25coU2bdpw/fp1QDcM9I8//sDd3T1HzzFuRfrkk08yXTtICCGEyG8kQSqApk+fzqVLlwDdNPfr16/H09Mzx89p2bKlOuHXkSNHWLdunVXjFEIIIWxFEqQC6ODBg+r2unXrMl3z53E0Go1JK9LUqVOfODYhhBDCHkiCVMAoisKJEycAKFmyJEFBQU/0vC5duhAaGgroppDfsWPHE8cohBBC2JokSAVMXFycWphtjansHRwcGDNmjLovrUhCCCGeBpIgFTDHjx9Xt6211k/Pnj0pXbo0AOvXr3/sQohCCCEKrn79+vHcc8/ZOozHkgSpgDF0r4H1EiRnZ2eT5UY++eQTqzxXCCGE7W3fvh2NRsPdu3dzdN+FCxfQaDTp/mieO3fuYxfMtQeSIBUwuZEgAbz22msULVoUgJUrVxIdHW21Z1uboiicO3cOrVZr61CEEKLA8fb2xsfHx9ZhPJYkSAVMbiVI7u7ujBw5EgCtVsv06dOt9mxrGzFiBOXLl6d79+62DkUIIdBqtUydOpWyZctSqFAhatSowS+//IKiKLRp04Z27dqp88zdvn2bkiVLMmHCBOC/1p21a9dSvXp13NzcaNCgAVFRUSavsXPnTpo2bUqhQoUIDg5m2LBhJCYmqucfPXrEmDFjCA4OxtXVlZCQEJYsWcKFCxdo2bIlAL6+vmg0Gvr16wfAhg0baNKkCT4+PhQpUoROnTpx9uxZ9Zlly5YFoFatWmg0Glq0aAGk72J79OgRw4YNIyAgADc3N5o0acL+/fvV84b3uHXrVurWrYu7uzuNGjXi1KlT1vkfkBlF5Mjnn3+uVKlSRalYsaICKPHx8bYOKUdKlCihAIqfn5+i1Wqt+uy7d+8qXl5eCqC4uLgoV65cserzrSEtLU3x9vZWAAVQYmJibB2SEKKA++ijj5TKlSsrGzZsUM6ePat8/fXXiqurq7J9+3bl8uXLiq+vrzJnzhxFURSle/fuSnh4uJKSkqIoiqJs27ZNAZQqVaoomzZtUo4cOaJ06tRJKVOmjJKcnKwoiqJER0crHh4eyuzZs5XTp08ru3btUmrVqqX069dPjeHFF19UgoODlV9//VU5e/assmXLFuWnn35SUlNTlVWrVimAcurUKSU2Nla5e/euoiiK8ssvvyirVq1Szpw5o0RGRiqdO3dWqlWrpqSlpSmKoij79u1TAGXLli1KbGyscuvWLUVRFKVv375K165d1dceNmyYEhQUpKxbt045duyY0rdvX8XX11e93vAe69evr2zfvl05duyY0rRpU6VRo0YWfb/j4+Oz9fktCZKFsvsNtieGmAGlcePGufIaY8eOVV/j7bffzpXXeBKnTp1S4wOUzz77zNYhCSFyS506ilKiRN5/1amT7RCTkpIUd3d35d9//zU53r9/f6Vnz56KoijKzz//rLi5uSljx45VPDw8lNOnT6vXGZKHn376ST1269YtpVChQsqKFSvUZ73++usmz9+xY4fi4OCgPHz4UP29uHnz5gxjNLzGnTt3snwvcXFxCqAcPXpUURRFOX/+vAIokZGRJtcZJ0j3799XnJ2dle+//149n5ycrAQFBSnTp083ef0tW7ao16xdu1YBlIcPH2YZU0ay+/n9+FVJxVPj5MmT6rY1u9eMjRgxgjlz5pCUlMTChQt599138fPzy5XXssS+fftM9levXs2QIUNsFI0QIldduwZXrtg6iixFR0fz4MED2rZta3I8OTmZWrVqAdC9e3dWr17NJ598woIFC6hQoUK65zRs2FDd9vPzo1KlSmpJxeHDhzly5Ajff/+9eo2iKGi1Ws6fP8/Ro0dxdHSkefPmOYr9zJkzTJgwgb1793Lz5k21rjMmJoawsLBsPePs2bOkpKTQuHFj9ZizszPh4eEmJSEA1atXV7eLFy8OwI0bNyhVqlSO4s4uSZAKkNyqPzJWrFgxXnvtNb744gsSExP5/PPP1b5ye2Dcrw3w999/c+vWLYoUKWKjiIQQuSYw0O5f9/79+wCsXbuWEiVKmJwzLB7+4MEDDh48iKOjI2fOnMlxOPfv3+eNN95g2LBh6c6VKlXK4kE1nTt3pnTp0ixevJigoCC0Wi1hYWEkJydb9LzHcXZ2Vrc1Gg1Arg62kQSpAMmLBAngnXfe4csvvyQtLY158+bx9ttv4+HhkWuvlxPmLUhpaWmsXbuWPn362CgiIUSuOXDA1hE8VmhoKK6ursTExGTagvP222/j4ODA+vXr6dChAx07dqRVq1Ym1+zZs0dtSblz5w6nT59Wf8/Xrl2b48ePExISkuHzq1Wrhlar5e+//6ZNmzbpzru4uAC635cGt27d4tSpUyxevJimTZsCukLwx91nrnz58ri4uLBr1y51Pr2UlBT279/PiBEjMr0vL8gotgIkrxKkMmXK0LNnT0D3Q7R48eJce62cSElJITIyEvjvBxd03WxCCGELnp6ejB49mpEjR7Js2TLOnj1LREQEn332GcuWLWPt2rUsXbqU77//nrZt2/LOO+/Qt29f7ty5Y/KcSZMmsXXrVqKioujXrx/+/v7qSLExY8bw77//MmTIEA4dOsSZM2f47bff1PKCMmXK0LdvX1577TXWrFnD+fPn2b59Oz///DMApUuXRqPR8OeffxIXF8f9+/fx9fWlSJEiLFq0iOjoaP766y9GjRplElNAQACFChViw4YNXL9+nfj4+HTv38PDg0GDBvHOO++wYcMGjh8/zsCBA3nw4AH9+/fPhe94DuS4ukkoipI/i7QrVKigAIq7u7s6yiC3REVFqYXQ5cuXz/XXy46DBw+qMb300ktKQECAAiiFChVSEhMTbR2eEKKA0mq1ypw5c5RKlSopzs7OStGiRZV27dop27dvV4oVK6ZMmTJFvTY5OVmpU6eO8uKLLyqK8l8B8x9//KFUrVpVcXFxUcLDw5XDhw+bvMa+ffuUtm3bKoULF1Y8PDyU6tWrKx9//LF6/uHDh8rIkSOV4sWLKy4uLkpISIiydOlS9fykSZOUwMBARaPRKH379lUURVE2b96sVKlSRXF1dVWqV6+ubN++XQGU1atXq/ctXrxYCQ4OVhwcHJTmzZsripJ+FNvDhw+VoUOHKv7+/oqrq6vSuHFjZd++fer5jIrEIyMjFUA5f/58jr/f2f381iiKfnIFkSMJCQl4e3sTHx+Pl5eXrcN5rEePHuHu7o5Wq6VWrVpERETk+mu2bduWLVu2ALr5Mtq1a5frr5mVhQsXMmjQIADmzJnDsWPH1Nat1atX54up74UQwtj27dtp2bIld+7cyReTL9qD7H5+SxdbAXHmzBm1mC03u9eMvfXWW+r2ggUL8uQ1s2JcoB0eHm6SEEk3mxBCCGOSIBUQeVV/ZKxz587qqIw//viDmJiYPHndzBgKtJ2cnKhZsyatW7emcOHCgC6+1NRUW4YnhBDCjkiCVEDYIkFycnLi9ddfB3RDMW1ZrH3//n2OHz8O6EZsFCpUCFdXVzp06ADoRn38888/NotPCCEs0aJFCxRFke61XCAJUgFhiwQJYMCAATg56WaTWLx4ca7Nj/E4ERERahdjvXr11OPdunVTt9esWZPXYQkhhLBTkiAVEIYEydHRMdO5MHJDUFCQWutz/fp1m9X6mNcfGXTo0EGdfGzNmjXImAUhhBAgCVKBkJaWpq56HBISYjIHUF6wh2Jt4wkijVuQvLy8aN26NQCXLl3Kk9F9Qggh7J8kSAXAxYsXSUpKAvK2e82gRYsWVK5cGdAt7XHs2LE8j8GQILm7uxMaGmpyzribTUazCSGEAEmQCgRb1R8ZaDQadf4h0M1HlJfi4uK4cOECAHXq1FFrogy6dOmirusjCZIQQgiQBKlAsHWCBNCnTx/c3d0BWLZsmbpAY14wrj8y7l4zCAwMVFfCPn78OKdPn86z2IQQQtgnSZAKAHtIkHx8fHj55ZcBuHfvHj/88EOevXZmBdrGjCeNlNFsQgghJEEqAIwTJEMtUI4oCsTHw5kzEB0NFy7AlStw/Trcvg337sHDh5Caqrs2E8bdbPPnz8+zEWOZFWgbK2gJUmRkJEeOHLF1GEIIYbecHn+JyM8URVETpODgYHXmaFVyMsTGwtWruqTH8GW8f/UqJCZm/0WdnMDZWfffQoWgWDEoXpzaxYvzTfHiHIyNJfbIEaIWLKDaM89A8eLg4WHFd/0fRVHUFqQiRYpQtmzZDK+rUKECVatW5dixY+zevZvY2FiKFy+eKzHZ2j///KNOLhcZGUnNmjVtHZIQQtgdSZCectevX+fu3buAvnvtxg347jv45Rc4e1a3b22pqbov0LUu3bgBR48C0Ff/BcDgwf/d4+mpS5SKF4fgYKhYESpV0n1VqAD6+qWcunjxInFxcYCu9chQjJ2Rbt26qSPsfv/9d9544w2LXtPeffnll2rr3ebNmyVBEkKIDEiC9JQ7ceIETkAH4KPoaChR4r/kJTs8PXX3BAVBYCBoNJCSonuG8X8zO5aYCNeu6bazcu+e7iuzAung4P8SJsNXxYpQqhQ4ZN5TbNy9lln9kUG3bt346KOPAN1otqcxQXr48CG///67um+YH0sIIYQpSZCeZlFR+H70EZeBYgDnzpmeL1lS9xUUpEuCDImQ8ban55PHoSi6WqXYWLh2jR9nzSJy/XqKA89Wr04VX1/dudhYXZKUkUuXdF9btpged3PTtTBVrQo1a/73VawYkL0CbYNatWpRqlQpYmJi+Ouvv4iPj8fb29vit22PNm7caDKCUBIkIYTImCRIT5s7d+DHH+Hrr+HAAWqanw8Kgj59oF8/XStMXtBooEgR3VdYGPXKlOHlChUA+OzePc5ERODo6Ki79v59XRH4qVPpv/RdhSaSknTdd0ePwk8//Xc8MBBq1KDm8eP0AA4B9WrXfkyYGp577jnmzZtHSkoK69ato2fPnlb4BljXiRMn2LNnDy+//DKurq45unflypUm+5IgCSFExjSKLD5lkYSEBLy9vYmPj8fLy8u2waSl6VpWvv4a1qyBR49MTj8CfgPa/fgj3i+8oCuetrFnn32WjRs3ArB27Vo6dOiQ9Q2KAjdv/pcsnT793/bZs4/vwgNdwXi1av+1MtWtC9Wrg1GSsW3bNlq1agVA9+7d+fnnny18h7kjMTGRkJAQrl27Rq9evVi+fHm273348CEBAQHp5qC6ffs2vr6+1g5VCCHsUnY/vyVBspDdJEg3bkCLFmA0lF9Vuzbjz55lYXw8miJFuHnzZp6Hl5nffvtNHVrfqVMn/vjjD8sflpKiS5QOHTL9unXr8fe6uOiSpfBwCA8ntXZtAps25dadOxQuXJi4uDjc3Nwsj83K/vzzTzp37gyAg4MDp0+fpnz58tm6d82aNSbLqhjs3r2bBg0aWDVOIYSwV9n9/JZ5kPK7DRtMkyN/fxgxAg4fJv6vv5gSH89tbDdBZGY6duxIcHAwoGtBMiwFYhFnZwgLg169YMYMXWtaXBwrZsygIzAeOFWtGoSEpL83ORn27YPPP4c+fXAKCyPm/n02AePu3yfqo4900xzYiQ0bNqjbWq2WWbNmZfte4+61tm3bqtvSzSaEEOlJgpRD8+fPJzQ0NNMJB/OccV3Oe+/p5i2aPRuqV+fkyZPqKXtLkJycnHj99dcB3VxFixYtsu4LaDRsj45mHTAFiJ03TzfRZUIC7NoFn32mq8XKYOJM95QU2gLvAnU//lhXsF6yJDz/vC4B+/dfXe2TDRgnSABff/11tloGjUev+fr68tZbb6nnJEESQoj0JEHKocGDB3P8+HGT0VE2ZTyBY926ui4jPXtYYiQrAwYMUBeO/eqrr3hkVjv1pAxD/DUaDXXq1NEd9PSERo1gyBBYtkzX+nbnjq7VacoUeO45tIGB6R925QqsXg3vvAONG4O3NzRsCG+/DatW6Ubg5bLo6GjOnj1rcuzhw4d88cUXj73XePRat27dCAsLU89JgiSEEOlJgpTfGRfcms1Gbe8JUmBgIM8//zwAcXFx/Prrr1Z7dlJSkrqURmhoKJ5ZTVfg4wOtW8O4cbB6NQ6xsQx49lm6AVOBO3XqgHk/dXIy7NkDs2bBCy/oRgeWLQuvvALz50NkZM7mm8qG9evXq9tvvvmmOvLv888/5+HDh1nea1xs3r17d8qUKYOzszMgCZIQQmREEqT8zrgFyWwZEXtPkACTrp7stIRk16FDh0jVJyiWdIc2f/ll1qDrZvuoeXNdK1NUFCxeDK++mvEUCRcuwA8/6Fqnatf+L/GaOFHXQmU2eiynjLvXBg0aRPfu3QFdcvntt99met/Dhw/VInhfX19at26Nk5MTIfqarOjoaNLS0p4oNiGEeNpIgpTfGSdImbQgubu7qwXR9qZZs2aEhoYCsHPnTo7qlyR5UjmZIDIjHTt2VFtoVq9ejaLR6CajHDAAli6FkychLg7++EPX8tS8uW4aAWOJifDXXzBpErRtq0uYwsNh1Chdd51+CZTsSEpKYtu2bQAUL16catWqMXr0aPX8zJkzM01yNmzYYNK9Zmg5qqRP8h49esTFixezHYsQQhQEkiDld5l0sSUlJXFOP3N25cqVcchiOQ5b0mg0DBo0SN1fsGCBVZ5rvMSIJS1Ifn5+tGjRAoDz589nnLj5+0OnTrrape3bIT4e9u+HuXOhRw9dYbextDTd+dmzdQXfAQG6IvGBA+Hbb3UznWcy68aOHTvUbrRnn31Wratq2bIlAGfOnDFZQsSY8eg1Q6sT/JcgAZzObIkXIYQooOzzU1NkXyYtSGfOnEGr1QL2271m0Lt3bzz0sS9fvpx7mS03kgOGFiQXFxeqV69u0TMM8zSBrhXpsZyddYXyw4bpZvW+dAkuXtQtDvzGG7oWKHOnTsFXX0HfvlC+vC6p6tsXfv/dZKSccffas88+q26/88476vaMGTPSPd589Frr1q3Vc8YJktQhCSGEKUmQ8rtMapDyQ/2Rgbe3N6+88goA9+/f57vvvnui5929e1f9wK9ZsyYuRiP7cqJr167q9po1aywLplQpXeH2woW6GqabN+G332D0aGjQIP2s5lev6lqTunaFokWhZ09YtYrta9cCuskh27Rpo17+7LPPUlWfeP3777/8+++/Jo/bsGEDifp/I8bdawAVK1ZUtyVBEkIIU5Ig5XfGCZJRDUx+SpAAk262L774gieZ4P3AgQPqtiX1RwbBwcHUrVsX0BV9nz9/3uJnqYoUgS5d4NNPYfduXbfctm3/1SkZ15Hdv69riXrhBXacOsVK4L3y5fEzSqo0Go1JLdKnn35q8nLGo9defPFFk3PSgiSEEJmTBCm/M9QgubuDUZ1RfkuQatasScOGDQGIiorin3/+sfhZT1qgbcx4aQ6LW5Gy4u6uWyrm/fdh0ybd8ih//qlbTNhofTR34AXgwzNndLVLXbrAihXw8CEvv/wyQUFBgG4JF0M9kfHoNT8/P3WNOQN/f3/8/PwASZCEEMKcJEj5naEFKZMh/sbDue2d8ZD/CRMmWNyK9KQF2sZyPUEy5+oKHTvqFh6+fh02bmRTmTLcML7m0SPd6LmXXoLAQFzeeovpnTqhQTcruWH5kfXr12favWZgaEW6cuVKukVshRCiIJMEKb8zJEhGXTNpaWlqi0BISEiGH4z2qEePHlSoUAGAf/75J92yGtllaEHy8vIyqbOxROXKldVn7Ny5k7hMhuZrtVoSExO5ceMG58+fJyoqir179xITE2P5izs7k9yiBS/cukUQ0NXLC+Wtt6B48f+uSUiAJUt4ZdEizms0TAJ2Ll3KjRs3Mh29ZkxGsgkhRMYkQcrvDH/1GyVIFy5cUJftyA/dawbOzs58/PHH6v64cePUkXjZdfXqVa5cuQJA3bp1n3h6A41Go7YiabVa2rRpQ/369QkLC6Ns2bIEBATg4eGBo6MjhQsXplixYpQrV45q1arRoEEDypcvz1p9gbUldu/ezb1790gDCnfqhGb+fLh8WTe/0quvmrQcllYU3geiUlJIqlGDwF9/xZuMu9cMpA5JCCEyJglSfqbVgmGJiXw6gs3c//3f/6nrph0+fJgff/wxR/dbs/7IwHi4/5EjR9i3bx/Hjh3jwoULxMXF8eDBg0zvTU1N5bXXXsu05elxjFvR2rdvr9twcICWLXUTVl6/Dt9/D88+i2KUDJa6do3ZyclcBX4vUgTnAwcynGNJEiQhhMiYJEj5mfEHs1ELUn5OkBwcHPjkk0/U/ffff5/k5ORs32/N+iOD8PBwOnXqZHLMw8ODgIAAypYtS1hYGPXr16dVq1Z07tyZl156iQEDBqivf+PGDV5//XWLaqqME6Rnnnkm/QXu7vDyy7B+PZrLl/m+dm0ijU8Djc+c0S3QW60azJunWzZFTxIkIYTIhCIsEh8frwBKfHy87YKIjVUUXbuAonTpoh5+9dVXFUABlAMHDtguvifQpk0b9T189tlnFt136dIlq8Z0584dJTExUdFqtdm6/tq1a0rRokXVeL7++uscvd7Vq1fVe+vUqZOtew4dOqQASjVQ5oJyV6P579+I4cvNTVF69VKUHTuUpIcPFQcHBwVQatWqlaP4hBAiP8ru57e0IOVnmcyibdyCVLly5byMyGqmTp2qbk+ePDlbI6y0Wq06B1Lx4sUpUaKEVWPy8fHB3d0djUaTreuLFSvGokWL1P1hw4Zx4cKFbL/exo0b1W3j2bOzUqNGDdq2bctRYDgwrk8fWLYMGjf+76KkJN3s3k2b4tqoEW8XKYIruiJt5TGtXEeOHOG7775Tlz0RQoinlSRI+VkGs2griqImSKVKlVKX8Mhv6tatq468unHjBrNnz37sPdHR0dy9exfQda9lN5HJTc899xz9+vUD4N69e/Tt2zfTRWXNZVh/lA0TJ07E1dUVV1dX+g8dCn36wM6dcOwYjBhhMr8SkZFMj4vjEjA2MZFYo0k2zd28eZOmTZvSu3dv2rZtS3x8fLZjEkKI/EYSpPwsgxaka9euqR9c+a3+yNxHH32Eo6MjoJsh+nGFzrlRoG0Nc+fOpXTp0oBu+oLsJHtpaWls2rQJ0C3FUr9+/Wy/XuPGjYmKiuLYsWNqwTsAoaG6hXINy5kY1WgVBd4DAhs21C1vsnt3uqLudevWkZCQAMCuXbto1aoVN2/ezHZcQgiRn0iClJ8ZdzvpE6T8XKBtrmLFigwYMADQtb5MmTIly+tzo0DbGry8vFi2bJnaojV+/HiOHj2a5T379+/njr6Yum3btjiZr9n2GCEhIZQvXz7jk25u0Ls37N0Lu3dzpm5dUvSnHNLSdMubNGqk65ZbvVo3WhJdgmQsIiKC5s2bExsbm6PYhBAiP5AEKT/LoIvtaUqQQDejdiH9GnNffPEFFy9ezPRa4xYkwxpq9qJ58+a8/fbbACQnJ9O7d291rqqMGHevZbf+KMc0GmjQgMvTp1MamAzcM1rPj9274fnnoUoV0hYsYLs+Ji8vL3Vpk+PHj9O0adMs/78IIUR+JAlSfpZBF9vTliAFBQUxfPhwQJdYTJw4McPrUlJSiIzUDXAPCQlR1xizJ5MnTyYsLAzQzfH0wQcfZHrt+vXr1e127drlalyVKlUiFpgAvNKkiW6ZE32cAJw+jeNbb3EoPp5xwAutW7Njxw7KlCkDwNmzZ2nSpInMxC2EeKpIgpSfPeVdbAZjxozBV19Y/O233xIVFZXumqioKJKSkgD7qj8y5ubmxnfffacu/TJ9+nR27tyZ7rqbN2+qrWHVqlWjZMmSuRpX8eLFKaxvgYyKjtYtlHvkCKxfr5uQUi8QmAIsXLeOcgsWsGvVKnUZlsuXL9OsWbPHdh0KIUR+ke0Eae/evbkZh7BEFl1s/v7++Pv72yIqq/Px8WHs2LGAbpTeu+++m+4a4/oje02QQDcMf9KkSYBuWoI+ffpw7949k2s2b96sDrfPte41IxqNRp0w8sKFC7pEU6OBZ5/VLWmyfz8bvb0xjL1zfvQIZswgqHFjIps1o41+Konr16/zzDPPqIXcQgiRn2U7QRo/frxJE/rWrVtp06ZNrgQlssmsiy0+Pl4tmH1aWo8Mhg4dqs5r9Mcff7Br1y6T8/ZaoJ2Rd955h8b6eYnOnz/PqFGjTM7nSf2RGUOCpCgK0dHRJucuBwbybHw8FYFfAgLA1VV3IikJ96++YtO5c6wqWpQy6EZRZmeUnhBC2LtsJ0jfffcdb7zxBt999x0NGjRg4sSJDBs2LDdjE49j1sX2NHavGRQqVMikZmfs2LEmkxoauqQcHR2pVatWXoeXI46Ojnz77bdqt9ZXX33FH3/8AehalQwJkoeHh5pI5baslhwx1EOdA4699RacPw+jRumWOQE0yck8HxfHGeAr4Ofp0y1ee04IIexFthOk6Ohopk2bxpgxY3j//ffZuXMnXbp0yc3YxOOYtSA9zQkSQL9+/dQP8p07d7J27VoAEhMTOXbsGKCr2SlkPBLLTpUrV86kpWXAgAHExcVx6NAhbty4AUCrVq1wNbTW5LKsEiTD9xmgQ4cOULw4zJwJFy7AuHHg6QmAE9AfiHzwgNNt2sDly3kQuRBC5I5sJ0ifffYZr732GvHx8bz77rv873//Y9myZbkZm3gcsxqkpz1BcnJy4uOPP1b3x40bR1paGhEREWj1c/XYc/2Ruf79+9O5c2dAN1v4G2+8YTJ6LSezZz+pzBKkR48esWXLFgACAgJMJ54sWhSmTNElSh98gNbLCwAXoPGRIyghITB8OFy7lhdvQQghrCrbCdKKFSuIiorizp07fPfdd9SqVUuG9dpaAWtBAnj++efVJCgqKooffvghX9UfGdNoNCxevFgtpl+9ejWffvqpej63h/cbq1ChgrptnCDt2LGDRP2/s/bt2+PgkMGvDD8/mDgRh4sX2dKwIYaSc82jRzBvHpQrB2PGgH4ZGCGEyA9yPIrN2dmZatWq0bNnT5O/5oUNZFKD5OHhQXBwsI2Cyl0ajYZPPvlE3Td09xrkpxYkSL+grWGZmIoVK1KuXLk8i8P438ypU6fU+i7j2bM7dOiQ9UN8fKj955/U8PRkGvDAcPzhQ5g+XZcozZoFWUyQKYQQ9sLiUWx//fWXjGKzNaMWpCQnJ86fPw9A5cqV7WKh1tzSsmVLnnnmGQAuXrzImjVrAHB3dyc0NNSGkVmmW7du6oK2Bnk1es2YoZvt7t27apG1of7I0dFR/Z5nxc/Pj/5jxjAWKAusrVDhv1Fvd+7A229D5crw/ffqEiZCCGGPLB7FNmHCBBnFZmuGBMnRkdMXLqh1OE9r95ox41Ykg9q1a+d4zTJ7YbygLdgmQTJM+ghw+vRpoqOj1T+KGjdujI+PT7aeM3z4cAICArgBdDpzhqO//KJb+82QtF+4AL16Qd26sHWrdd+EEEJYiYxiy88MXWweHpw4eVI9XBASpFq1avHSSy+ZHMtv3WvGvLy8+OWXXwgNDaVHjx7Zaq2xNvNCbeOC8cd2rxkpXLgw77//vro/+rPP4NtvISICjN9XZCS0aaObkPLwYXbs2EHDhg3p378/qampT/ZmhBDiCVllFNuIESP4+uuvOXjwYJYLcAorM7QgFYARbBmZPHmySYtRfirQzkjdunU5duwYP/30E46Ojnn++uYJUo7qj8y8/vrr6lptmzZtYvv27VCzJmzcCJs26bYNNm5EqVWLC82bc2XPHpYuXcrPP/9s+RsRQggrsMootlatWnHlyhWmTZtG7dq1CQsLo0ePHnz88cfqBHgiFxgSpAIygs1cSEgIgwYNAnRFxs2bN7dxRPmbcYIUGRnJtm3bAAgODlYX2c0uFxcXPvzwQ3V/8ODBbN26VVf83bYtHDwIy5ejlCoFgEZR6K0onAY+ARZMm2YyEagQQuQ1jZILv4WSkpKIioriyJEjHD169KlceiAhIQFvb2/i4+Px0s//kudcXCAlBWrWpHpaGkePHsXJyYkHDx6oC6I+7dLS0vj222+pUqUKDRo0sHU4+ZpWq8XDw0Nd9NfgjTfeYOHChTl+XlpaGjVq1FAn8QSoWrUqQ4cOpVevXiQnJ9PrhReo/NdfvAf4Gt0bC9wZM4bQKVMgo6kFhBDCQtn9/M5RgnT8+HF+/PFH3n777WwXbD6tbJ4gpaToEiRAadyYQgcO8OjRI6pUqcLx48fzPh7xVKhevTpHjx41Ofbbb79ZXG8YGRlJ165duXTpkslxHx8fvL29uXjxIgABTk5sbt2aqlu34mhcf9SggW4upVzuPlUUhSlTpnD16lWmTZumLgMjhHj6ZPfzO0d/mk2dOpWoqKgMk6OkpCROGhUKi1xmNMT/gYODWvtVULrXRO4w7mYDXVdZ69atLX5erVq1OHv2LCtWrDBZV+7u3btqcuTv788vf/1F9Q0bUKKi2GS8VMyePVC/PgwYAPolWHLD9u3bee+99/jiiy9M5qUSQhRcOUqQ9uzZk+nQfjc3NwYOHMjUqVOtEph4DKMEKcHoL25JkMSTME+QWrRogYeHxxM909nZmRdffJGdO3dy4MAB+vbti4u+9bN69ers37+fpk2bAuBUqRInpk7lGUCtqlMUWLIEKlSA2bN1radW9u+//6rbBw8etPrzhRD5T44SpMuXLxMSEpLp+TfffJPff//9iYMS2WA0i/Yto5GDkiCJJ2GeIOV09Nrj1KlTh2+++YZLly6xZcsW9u3bp452M3jttdfY7+NDdWC0gwNa/WK4JCTAqFFQowZs3mzVuCIiItRt6aIWQkAOEyQ/Pz9iY2MzPR8eHk50dPQTByWywagF6YbRtiRI4knkdoJkEBAQQOvWrXE1zLJtxNPTkzfffJNUYKZWyyevvgr9+/830eSJE7r5lLp1g5gYq8RjnCCdPHmStLQ0qzxXCJF/5ShBatasGd98803mD3NwSDcCxp7Nnj2bqlWrEhoayrBhw/LXsGKjpOhqQoK6bf4BJ0ROGP/7qVChgskitnlp6NCh6kjMT7/9lvtz5sDevbqibYM1ayA0FGbOfKJut9u3b3PhwgV1PykpyWRfCFEw5ShBGj16NIsXL860iHH37t15usDmk4iLi+Pzzz/n4MGDHD16lIMHD7Jnzx5bh5V9Rl1sl27dAqB06dJPXC8iCjZvb2+GDx+Or6+vyTxGeS0oKIhXXnkF0BV0L126VDeSbdcuWLYMihXTXZiYCKNH65Yt2b3botc6dOhQumPGUxMIIQqmHCVIderU4YsvvuCtt96ibdu2rFmzhpiYGG7fvs1vv/3GmDFjePnll3MrVqtLTU0lKSmJlJQUUlJSCAgIsHVI2WfUgnQrORmQ7jVhHXPmzOHWrVv07NnTpnG8/fbb6vbs2bN1y484OECfPnDyJAwe/F+325Ej0KgRvP463L6do9cx7l4zkDokIUSOZ2AbOHAg27dvJz4+nueff56yZctStGhRunXrRrVq1Rg5cqRVAvvnn3/o3LkzQUFBaDQadcV2Y/Pnz6dMmTK4ublRv3599u3bl+3nFy1alNGjR1OqVCmCgoJo06YN5cuXt0rsecIoQTJsSYIkrEVjSDxsKCwsTF2098KFC2w1XtjWxwc+/1zX7Var1n/HFy+GypVJ+Pxz2j3zDP37939st78kSEKIjFg0RW2TJk3Yt28fx48fZ/ny5Xz55Zf8+++/bNiwQR2++6QSExOpUaMG8+fPz/D8ihUrGDVqFBMnTiQiIoIaNWrQrl07bhjNlVKzZk3CwsLSfV29epU7d+7w559/cuHCBa5cucK///7LP//8Y5XY84RRF5shQbJVvYgQuaVXr17qdmRkZPoL6tWDfftg7lwwjHaLi8Nr6FDGbt7Mv0uX8vHHH2f5GoYEycXFBQf9rN2SIAkhnB5/SeYqV65M5cqVrRWLifbt29O+fftMz8+aNYuBAwfy6quvArBw4ULWrl3L0qVLGTt2LJBxbYHBypUrCQkJwc/PD4COHTuyZ88emjVrluH1jx49MlmIN8GoMNomjFqQDKmS4b0I8bQwXgMu07ogJycYNgz+7/9g5EhYuRKAlsAhYOqUKRz/v/8j1HiBXL179+6pa0rWqFGDO3fuEB0dzYkTJ9BqtWrCJIQoePLlT39ycjIHDx6kTZs26jEHBwfatGnD7mwWagYHB/Pvv/+SlJREWloa27dvz3IE2NSpU/H29la/goODn/h9PJEMutg8DX9BC/GUqFSpkpqkREVFZX1xiRLw889o//yTC/ouQlfgA60Wt0aN0BpNBmlw+PBhdfRq7dq1CQ0NBeDBgwfEWGkKASFE/pQvE6SbN2+SlpZGMcNIFr1ixYpx7dq1bD2jQYMGdOjQgVq1alG9enXKly+f5XpT48aNIz4+Xv0yX1sqz2WQINls0Vwhcombm5s6OW125yc6W7EiVRSFKYBhjvlyDx+iadIEhg836Z42rj8yTpBAutmEKOjyZYJkLR9//DEnTpzg2LFjzJs3L8vCVFdXV7y8vEy+bCqDGiRpQRJPI0M3W1JSEufOnXvs9QcPHiQJGA8Mrl+fA/rjGkXRLXxbtSps2ABIgiSEyFy+TJD8/f1xdHTk+vXrJsevX79OYGCgjaLKYxnUINk8aRMiF1StWlXdzs78RAcOHFC3u7z/Pgv69GE08MBwMCYG2reHXr04px/56ujoSFhYmCRIQgiVxUXaSUlJHDlyhBs3bqDVak3OZdVVZQ0uLi7UqVOHrVu38txzzwGg1WrZunUrQ4YMydXXthtSgyQKCPMEyfAznxnjxWbr1KlDgwYNqLxuHatv3mQR0Npw8vvvWQUMAU6GhuLm5mYy6EQSJCEKNosSpA0bNtCnTx9u3ryZ7pxGo7HKOkb37983Wdft/PnzHDp0CD8/P0qVKsWoUaPo27cvdevWJTw8nDlz5pCYmKiOanvqSRebKCCME6THFWprtVq126xEiRJqi/Ls2bPp3bs3bYAh7u7McXLCMSGBosAKYH9CAty4gUdAAGXLluX8+fMcP34cRVHU32kxMTGUKVPGLuaIEkLkAcUCISEhyltvvaVcu3bNktuzZdu2bQqQ7qtv377qNZ999plSqlQpxcXFRQkPD1f27NmTa/GYi4+PVwAlPj4+z17TRKNGigKKAoojKK6urraJQ4hc9ujRI8XJyUkBlGrVqmV57enTp9XfFV26dFGPa7Va5f/+7//Ucx1q11bO1Kih/gwpoCj+/oqyYoXSsWNH9bpLly4pWq1W6dGjhwIob775Zpavf/nyZWX8+PHKzp07rfHWhRC5ILuf3xYlSJ6enkp0dLRFgeV3n3/+uVKlShWlYsWKtk2Q9L/ckzQaBVD8/f1tE4cQeSA0NFQBFBcXFyUlJSXT63744Qc1ufnwww9Nzt25c0cpU6aMet7Ly0vpDkqccZIEyuGKFRV//TUbN25U1q5dq97j6+ubZZzt27dXAMXPz09JSkqyynsXQlhXdhMki4q0X3jhBbZv327Jrfne4MGDOX78OPv377dtIPouNhniLwoCQzdbcnKySde7OeP6o7p165qc8/HxYcWKFTg7OwO6yV5XAlWB1K5d1euqnz7NceAFdPMkjRo1Sj13584d4uPjM3zt2NhYNm7cCMDt27c5ceJEDt6hEMLeWFSD9Pnnn9O9e3d27NhBtWrV1F84BsOGDbNKcCIL+iLte/pJ7qT+SDzNqlatykr9DNnHjh3LdAZ/8wJtc+Hh4UybNs0k6fGrXBmn1athxQoYMgRu3aIosBJY/d576mLQBufPn6dmBrNy//zzzyYDVg4fPpzhdUKI/MGiBOnHH39k06ZNuLm5sX37dpOiRY1GIwlSXtAnSNKCJAoC80Lt//u//0t3jXmBtvlEsgYjRoxg27Zt/PHHHwDUqlULNBp46SVo2ZKUgQNx1p/rlpxMY2AQ8Kv+/swSpB9//NFk//Dhwzl7k0IIu2JRF9v48eP58MMPiY+P58KFC5w/f179ys5EbuIJKUq6BElakMTTLDtzIUVHR6trJGbUemSg0Wj45ptvCA8Pp3DhwgwdOvS/k8WK4fzbbwzx8+OW/lAAsAr4ESgCGf6OO3fuHHv37jU5ltVakMauX79Oamrq4y8UQuQpixKk5ORkevToIQs52kpSEuib8qUFSRQEISEhuLi4AJknSFnVH5nz8/Nj9+7d3L17l4YNG5qe1GiIrlePqsBqo8MvAVFAoW3b0j3vp59+SnfMeJ23zPz0008UL16cunXrSpIkhJ2xKMPp27cvK1assHYsIrsymEVbWpDE08zZ2VldTPr06dMkm9UFwePrj8w5ODjg6OiY4bmwsDCuA88Dq194Aa2PDwCBwFtr18Kbb5rMRWbcvWaYjfv27dtcuXIlyxhWrlyJoigcPnz48YvxCiHylEU1SGlpaUyfPp2NGzdSvXr1dEXas2bNskpwIhOyUK0ogKpWrcrRo0dJTU3lzJkzJt1ukPMEKStvvvkm27Zto0yZMjz77bdo7txhfXAw7Q1F2F9+CVu3wvLlRBUurCY3DRs2pFmzZuos3IcOHaJkyZKZvs7t27fV7cxqm4QQtmFRgnT06FFdYSOPn9n2aTN//nzmz59vldnCLSazaIsCyLwOyXhfq9WqCVJWBdrZFRISYpJwUagQoypUYM2pU8wCPACio6FxY+IaNsQJSAV69uxJ0aJF1dsOHz5Mp06dMn2dO3fuqNtSvymEfbEoQdqWQR98QTF48GAGDx5MQkIC3t7etglCFqoVBZD5SLYXX3xR3Y+OjubevXvA4+uPLFW2XDkWnTrFX8CxOnVwOXgQtFpa7trFHqCPRsOLL75o0ir0uEJt8xYkIYT9sChBmjRpUqbnNBoN77//vsUBiWyQhWpFAZTVSDZrdq9lpmzZsgBEAwdnz6bh33+jfPghmtRU6gARGg2uK1ZQ5M03cXNzIykpSR3qrygKf/zxB4GBgYSHh6vPNE6QpAVJCPtiUYK0evVqk/2UlBTOnz+Pk5MT5cuXlwQpt2XQxSYtSOJpV758eVxdXXn06JFNEySA85cu0fC99/i7UCGKjR5NFcBVq4Xhw3H6/XdaV6rE2sOHiY6O5v79+3z11VeMHDkSBwcHzp49S5kyZUhOTibR6I8dSZCEsC8WJUiRkZHpjiUkJNCvXz+6dev2xEGJx5AuNlEAOTo6UqVKFQ4dOkR0dDSPHj3C1dUVgAMHDqjX5VaCVK5cOXXb0B32V3w8nwKfAMMNJ7duZYWrK32BVYrCrl27mDx5MqCrlTp06BBlypQxqT8yPFOr1cr0KULYCav9JHp5efHhhx9K61FekC42UUAZutnS0tI4deoUALt27WLXrl2AdQq0M2PSgqRPkA4ePEgSMAKI/fZbKFECAI9Hj/gFWAy8M2iQSVdaXFwcYNq9Brr55WJjY3MldiFEzln1T5X4+PhMF3IUViRdbKKAMi/Uvn79Oi+++KI6yeLAgQNz7bUzSpAMS5v4+voS2KsXHD0KRsXjA4Bfzp+nttFzbty4AaRPkEC62YSwJxZ1sc2bN89kX1EUYmNjWb58Oe3bt7dKYCIL0oIkCijjBOnw4cMsXryYq1evAtCyZUvGjx+fa6/t4+ODj48Pd+/e5fz581y9epVr164BULt2bd2alL6+8NNPPGzenLTBgykMVAR2A+8BM8i8BQl0CVLTpk1z7T0IIbLPogRp9uzZJvsODg4ULVqUvn37Mm7cOKsEJrIgNUiigAoLC1O3586dy6NHjwAICgrixx9/xMnJol9p2Va2bFkiIyOJiYlh37596nGTuieNhkJvvUXrTz7hk0uXqAe4ANOBZ4A1Fy8CpKtBgv9akNLS0rh9+7bJnEpZuXr1Kh07dqRIkSKsXbtWrc0SQljOot8mBXm+DruYKFJakEQBVaZMGdzd3Xnw4IGaHDk5ObFy5cpcqz0yZkiQ0tLS+O2339TjGRWGe9etS6NLl/gQGIuunqENEL52LaxZk2EL0vnz50lNTaVhw4YcPHiQZcuW0bt378fGtXz5cnXOpQ0bNtC1a1fL3qAQQpXjGqSUlBRat27NmTNnciMeuzd48GCOHz/O/v37bReEWQ2Sm5tbuuVehHgaOTg4UKVKFZNjs2bNolGjRnny+sZ1SMYJUu3atdNd2717d1KBz4sX5+6qVVzWH/dKSYFu3Wi8fDnuZvecO3eOnTt3cuDAARRF4bvvvstWXIZuRuCx678JIbInxwmSs7MzR44cyY1YRHaZdbFJ65EoSKpVq6Zu9+zZkyFDhuTZaxsnSIYuMm9vb8qXL5/u2p49e3Ly5EmioqLwe/552gUG8ovR+XoRERwAahodO3fuHH/++ae6n93WekPhN6DWRQkhnoxFo9h69erFkiVLrB2LyC6zLjapPxIFyciRIylbtiydO3dm8eLFuuLoPGKcIBnUqlUr0xgqVaqEn58fAM7FitEdeN3REcVd13ZUBdgLvO/hgQaIjY1l1apV6v0XLlzIVne+JEhCWJ9FNUipqaksXbqULVu2UKdOHTw8PEzOz5o1yyrBiUyYdbGVkBYkUYBUr17dZsPhM0qQsjsxZUBAAACL09KYsX0719u0oUJCAi7ApMREwoF+6JIig5SUFK5evUpwcHCWz75+/bq6LQmSENZhUQtSVFQUtWvXxtPTk9OnTxMZGal+PW5xRmEF+hYkLfAQaUESIq+UKVMm3bHsJkjGI9Ku+/jQr2JFphmd7wQcAswH+RuSwa1bt9K+fXt+//33dM82bkEyTpaEEJazqAVp27Zt1o5D5IQ+QXoAKEgNkhB5pVChQgQGBpq00mRUoJ0R4wQpLi6OG3fvMhbYV7gw3zk4UCghgZLANuADYAq6P4LOnTtH8+bNGTJkCCdPnuTIkSN06dJFfVZaWho3b95U96UFSQjrsKgFKSYmBkVRMj0ncpm+i01m0RYi7xl3s3l6elKhQoVs3WecIN24cUMt8j4UEMDhb75hq/6cIzAZ2AQEoivUvn//PidPngR0I9aM51C6deuWye/ja9euZfr7WQiRfRYlSGXLllVngzV269atDPvohZXpW5Bkkkgh8p7xorW1atXK9uKyhhokME2Q/Pz8KFmvHs8A7wOGkuzWwGHAY+dOjh8/bvKs6Ohok2cZS05OliWfhLACixIkRVEyHLVx//593Nzcnjgo8Rj6BEkmiRQi7xn/EZjd7jUwbUGKjo5Gq9UCugQpKCgIJxcXPgJaAg/1I98CgDHbtuH6wQcm9RBZJUgg3WxCWEOOapBGjRoFgEaj4f3338fd/b9pztLS0ti7dy81a9a0aoD2xuYzaaelwcOHgHSxCWELxr/jmjVrlu37jBOkU6dOqdt+fn44ODhQsWJFoqKiOOjuDocOsblCBdrqZwuvsX49fwM9gRiylyBVrlw527EJIdLLUYIUGRkJ6FqQjh49iouLi3rOxcWFGjVqMHr0aOtGaGcGDx7M4MGDSUhIwNvbO+8DePBA3ZQWJCHy3nPPPcdHH32Ek5NTjpb0MO5iM0+QAKZNm8ZHH33EoEGDKBQczIe1a7N+926mAc5AI3Sj3F5FWpCEyAs5SpAMo9deffVV5s2bJx/MtiAL1QphU46OjowfPz7H9xm3IJ09e1bd9vX1BaBDhw506NBBPV62fHlm797NTuAnoBzgC6wBVm7aBI8egatrhglSdob6JycnM3XqVI4ePcqMGTMynMJAiILMohqkChUqsHLlynTHly5dyrRp0zK4Q1iNLFQrRL7k7e2trpmYmpqqHje0IJkz1DrtB2oDxr9xu1+7Bo0awdmzGSZDj2tBunz5Ms2aNeODDz5g1apVTJo0KUfvRYiCwKIEadGiRRn2b1etWpWFCxc+cVAiC2azaIO0IAmRH2g0Gvz9/dMdzyxBMh4tFw+8CAwCkgwHIyKgTh3KHD6c7t6sEqTExEQaN27M3r171WM7duzIxjsQomCxKEG6du0axYsXT3e8aNGixMbGPnFQIgsZdLFJC5IQ+YNxHZLB41qQjH3j5kZ9QK1gio9n/P79TEM3f5JBVgnSpk2b0s1XFx0dbTLZpBDCwgQpODiYXbt2pTu+a9cugoKCnjgokYUMutikBUmI/MG4DskgOy1IBs888wxHgHrApQYN1OP/A/52dKSkoy5NyqoG6ciRI+p26dKl1e09e/Y8JnohChaLEqSBAwcyYsQIvv76ay5evMjFixdZunQpI0eOZODAgdaOURjLoItNWpCEyB8ySpAMRdrmgoKCTEYKA+qouXvA8k6dYM4cUvTnGqelcVBRaE7WLUhHjx5Vt19//XV1WxIkIUxZlCC988479O/fn7feeoty5cpRrlw5hg4dyrBhwxg3bpy1YxTGZBSbEPlWTlqQHB0dTVp4ypQpQ61atdT96LNneTBwIM2AS/pjAVotW4F+166hNSoEN2ZoQSpUqBC9evVSj+/evTtnb0aIp5xFCZJGo2HatGnExcWxZ88eDh8+zO3bt5kwYYK14xPmZBSbEPlWRjVImbUggWk3W1hYGCEhIep+dHS07ncwulFuh/XPdgSmKAopnTqB0ZptAA8ePFDnUKpatSqlSpWiZMmSAOzbt892E+AKYYcsSpAMChcuTL169QgLC8PV1dVaMYmsmCVIhQoVwskpR9NZCSFsxLwFyd3dPcvlmYwLtatVq4anpyfFihUDdAmSYQ6km8Di55/nt5o10eqvd924EWrXhoMH1WccO3ZMXci2WrVqADTQ1zLdv3+fY8eOPdH7E+JpYnGCtGPHDnr16kWjRo24cuUKAMuXL2fnzp1WC05kwKwGSbrXhMg/zBOkzLrXDCpVqqRuG5Y4MbQixcbGcu7cOfV8QPHi7OvQgfboEiYALlzQzZe0aBHoV0AwqF69OvBfggTp65BWr17NhAkT1IV1hShILEqQVq1aRbt27ShUqBARERE80q8XFB8fz5QpU6waoDBjVoMk3WtC5B/mCVJW3WsAffv2pX379rz00kt069YNwKSbzbhuKCAggGLFirEJXZdbXPnyuhPJyfDGG9CvHyeMWpMMLUgNGzZUjxknSKdOneL5559n8uTJzJgxI0fvU4ingUUJ0kcffcTChQtZvHixOjMsQOPGjYmIiLBacPZo/vz5hIaGUq9ePdsEYNbFJi1IQuQf5jVIj2tB8vX1Zd26dfz444/q71rjBMl4upWAgAACAwMBXdH28gEDYNiw/x727be8/u23GMq+DS1IVatWVS8xnh9pyZIl6rb84SsKIosSpFOnTmW4irW3tzd379590pjs2uDBgzl+/Dj79++3TQBmXWzSgiRE/pHTLraMGHe7HThwQN02TpAArt68CXPnwk8/gYcHABXu3+cA8IKPjxqLl5eXOp2A8bpuxs+WGlNREFmUIAUGBpqsJm2wc+fODCc3E1Zk1sUmLUhC5B8+Pj4mgyosSZDat29PcHBwuuPmCdK2bdtISUmBHj24tXYtD/Wj1fyBn+7ehRkzQFHQaDRqshQXFwfAvXv3TOpJCxcunOM4hcjvLJ4ocvjw4ezduxeNRsPVq1f5/vvvGT16NIMGDbJ2jMKYWRebtCAJkX8YJyPw+BqkjBQuXJgvv/wy3fGAgABCQkLUdTIjIiL44IMPePvttynasiXFL19mrf5aR4B33oGePSExUe36u3nzJlqtli1btuiSK71bt26RlKRbBS41NZU1a9aYzMgtxNPIovHhY8eORavV0rp1ax48eECzZs1wdXVl9OjRDB061NoxCmMyik2IfM14zUpLWpBA14rUt29fli1bBoCLiwve3t5oNBqWLVtGo0aNSEtLM6kdigc6Ax8A6ox1K1bA8eNU9/YmEl3yc/fuXdatW5fuNa9cuUL58uX58ssvGTJkCO7u7pw5c0aWlxJPLYsnihw/fjy3b98mKiqKPXv2EBcXx+TJk60dnzCnb0FKAZKRFiQh8hvjFiRLEySA2bNnq11qlSpVQqPRABAeHs4HH3xgcq2TkxM9e/ak2/PPo504Ee2qVWD43XH0KPP37aOd/tobN25kmCBduqSbr3vTpk2AbtLJ9evXWxy/EPbuiWYYdHFxoUqVKgDqD6fIZfoESRaqFSJ/slaC5Ovry9atW1myZAl9+vQxOTdu3Dj++usvtm3bRrFixfjll19o0qSJ6QNCQ+G55+DUKTySk1kHvAccPHCAq1evpnu9y5cvA5hMJrlt2zb69+9v8XsQwp5ZPFHkkiVLCAsLw83NDTc3N8LCwvjqq6+sGZvIiFmCJC1IQuQvxkP9nyRBAggNDWXmzJnUqFHD5LijoyObNm1i8+bNnDp1Kn1yBFC5MuzbB/oFcB2AKUC1Dz/EUJJdokQJ9fLLly+TmJhoMjnlX3/9pc7MLcTTxqIEacKECQwfPpzOnTuzcuVKVq5cSefOnRk5cqSsx5bb9DVI0oIkRP7UvXt3XF1dKVu2LI0aNcq113FycqJNmzZ4e3tnfpGXF/z6Kwe6dlWXKKkeHc2/QBlMJ5G8fPkyJ06cMEmIYmNjOX36dG6EL4TNWdTFtmDBAhYvXkzPnj3VY126dKF69eoMHTqUSZMmWS1AYUbfgmQo1ZYWJCHylyZNmnDt2jU8PDxMJtq1GQcHrr72GhN/+43vAR+gGrAf+CsggF/0l12+fJmoqKh0t//1118mczMJ8bSwqAUpJSWFunXrpjtep04dUlNTnzgokYnkZNAPvZUWJCHyLx8fH/tIjvSKFi3KOiAcOKU/5g+88OWXvKmvL80sQdq2bVtehSlEnrIoQerduzcLFixId3zRokW88sorTxyUyITZHEggCZIQ4skZ6qLOAPWBDfrjDmlpLFAUPgNiY2JMEiTD7Nrbtm2TOiTxVLJ4FNuSJUvYtGmTuhL03r17iYmJoU+fPowaNUq9btasWU8epdAxm0UbpItNCPHkjEfWxQOdgE+BkfpjQ4AqcXEM0i8/4uvrS/Xq1fn777+5efMmiYmJMtu2eOpYlCBFRUVRu3ZtAM6ePQuAv78//v7+Jn9hyNB/K5MWJCFELvD09MTFxYXk5GQA0oCZJUowctIkUgcMwElRaA2su3WLLoB/WJjJaLybN29KgiSeOhYlSNLnbCNms2iDtCAJIZ6cRqMhICBAnesIoGTJkvDaa3yxcSM9fv6ZYkAIsAdY7unJMX9/9dqbN29SpkwZABT9+m5C5HcW1SA9fPiQBw8eqPsXL15kzpw56gyrT7P58+cTGhpKvXr18v7FpYtNCJFLjLvZAHVB3NT69akHROqPewFvrV9PV6Ph/YZFbiMjIyldujSdO3dGq9UiRH5mUYLUtWtXvv32WwDu3r1LeHg4M2fOpGvXrhkWbz9NBg8ezPHjx9m/f3/ev7hZF5u7uzuOjo55H4cQ4qljniCVLFlS/e8loAmoQ/41ikK7rVtZjK4b4ubNm4BuoM6lS5f4888/OaCvVxIiv7IoQYqIiKBp06YA/PLLLwQGBnLx4kW+/fZb5s2bZ9UAhRGzBEnqj4QQ1mJcUwT/tSA1bNgQZ2dnHgBfP/ssiaNHq9cMANYDCTExAMTo/wtkuFyJEPmJRTVIDx48ULt2Nm3axPPPP4+DgwMNGjTg4sWLVg1QGDGqQbqPdK8JIawnsxak4OBgIiMjuX79Oi1atMDBwQFq1yatb18cU1JoA1SfPRteftmkhun69et5Gb4QVmdRC1JISAhr1qzh0qVLbNy4kWeeeQbQrQItrRq5SFqQhBC5JLMaJICqVavSqlUrXXIE0LMnZxctIk5/PuDWLahfn8Dz59V7rl27ltshC5GrLF6LbfTo0ZQpU4b69eur6/Vs2rSJWrVqWTVAYcQsQZIWJCGEtWTWxZYZj7ZtaQCcMByIi+O3e/d4Ub8rLUgiv7Ooi+2FF16gSZMmxMbGmqwi3bp1a7p162a14IQZs2H+0oIkhLAW4xYkBwcHAgMDs7y+SJEinAMaAVu8vakTH48bsAIoD5zStyBptdr/Wp6EyEcs/lcbGBhIrVq1TP7hh4eHU7lyZasEJjJgNsxfWpCEENZinCAFBQXh5JT1389ubm4ULlyYu0C/YsWIbd9ePTcF6PfPPzRr0IASJUpkuIabEPZO0vr8RGqQhBC5xLiLzVCg/TiGpCr21i229uzJWKNznW/d4sO9e3l07RrLli2zZqhC5AlJkPITsy42aUESQlhLyZIlKVasGIBaV/o4/vrZtG/fvs3FmBimAd2Bh/rzLYHdQPKJExk/QAg7ZvFitcIGzLrYpAVJCGEtrq6ubN26lT179vDiiy8+/gb+a0FSFIUjR44AuskkY4DfgWJAJWDipk2wfz/YYgUCISyU4xaklJQUWrduzZkzZ3IjHpEV6WITQuSiqlWr0r9//2y3Tvsbrcd26NAhdXsfUB84pt/3S0mBFi1g3boMn5OYmMhLL71Enz591AVzhbC1HCdIzs7O6l8KIo9JF5sQwo4YF3afNlqbDeAiuuVJ/jYcePAAunSBJUvSPWfy5MmsWLGC5cuXs3Tp0lyLV4icsKgGqVevXizJ4B+5yGX6FqSHgBZpQRJC2JZxC1JG7gLtgJ8NB9LSYMAA+PBDUBT1OuPPk40bN1o7TCEsYlENUmpqKkuXLmXLli3UqVMHDw8Pk/OzZs2ySnDCjD5BMnS0SQuSEMKWHpcgATwCXgJa9eqF/3ff6Q5+8AFcvgwLFoCTEwkJCer13t7euRKrEDllUYIUFRVF7dq1gfTNqhqN5smjEhkzS5CkBUkIYUvmy5MAuLu78+DBA5NjChDRuzfP1K4No0bpDn71FVy9CitWmNQdpaWl5WbIQmSbRQnStm3brB2HyA59DZK0IAkh7EFGLUg1a9bk33//TXf8ypUrMHIklCgBvXtDcjKsW4e2eXOKgrqumyxRIuyFxfMg7dixg169etGoUSPdP3xg+fLl7Ny502rBCSOKorYgGUq1pQVJCGFLGbUgderUKcNrr169qtt48UXYtAn0XWkOERHsBkL018kit8JeWJQgrVq1inbt2lGoUCEiIiJ49OgRAPHx8UyZMsWqAQq9pCS1qFFakIQQ9sC8Balz586MGDEiw2sNf0gD0Lw57NoF+hm7ywP/AuFIC5KwHxYlSB999BELFy5k8eLFODs7q8cbN25MRESE1YKzR/Pnzyc0NJR6eT3hmdkQf4DChQvnbQxCCGHEx8eHcuXKAbo5lL7//nsKFSqU4bVqCxK6OZNGf/01p5Yt44Z+9u6iwDag3o0bpKam8u677xIeHs7Bgwdz+20IkSGLEqRTp07RrFmzdMe9vb25e/fuk8Zk1wYPHszx48fZv39/3r6w2SzaHh4eODo65m0MQghhxMHBgbVr1zJ//nx27typtmobWpFGjhypLmhuaEE6e/YsTZs2ZebMmfR9910mtmrFVv3z3IE1wLGxY5k6dSr79+9nlKGoW4g8ZlGCFBgYSHR0dLrjO3fuVP+aEFYms2gLIexQ5cqVeeutt/Dx8VGPzZ49m1u3bjFr1iwCAwMBXYKUnJxMz549ua9vEd+7dy/7T5+mPfCT/l4noMbMmbyt3//nn39QjOZMEiKvWJQgDRw4kOHDh7N37140Gg1Xr17l+++/Z/To0QwaNMjaMQqQWbSFEPmKn58fAEFBQQDExsbSo0ePdK3vBw8eJAV4GfjM6PgMYJp+23w6GSHygkXD/MeOHYtWq6V169Y8ePCAZs2a4erqyujRoxk6dKi1YxQgC9UKIfKlEiVKcODAAQDWrFmT6XUKMAy4DnykP/Y/IADY/c8/VKpUSXedorB48WLu3bvHiBEjpNRA5BqLWpA0Gg3jx4/n9u3bREVFsWfPHuLi4pg8ebK14xMG0sUmhMiHGjdubLLv6urKwIED010XGhoKwMfA64Bhush+QO3Jk3VruaEr5XjjjTcYPXo08+fPz7W4hbB4HiQAFxcXQkNDCQ8PlxFVuc0sQZIuNiFEfvD222/z888/06tXL5o3b866desynAqgffv26vZioDu6ZUoAql+6BG3bwu3bbNq0Sb1u+PDhUp8kco1FCVKfPn1YunQpZ8+etXY8IjNmNUjSgiSEyA8cHBzo3r07y5cvZ/v27bRq1YqQkBCcnEwrPDp06GCyn9C6NSMqVUJdpe3ff6FZMzTG8ykB+/bty8XoRUFmUYLk4uLCJ598QoUKFQgODqZXr1589dVXnDlzxtrxCQOzGiRpQRJC5FcuLi6EhISo+76+vtStW9fkmiFDhuDQujXNAXVu7WPHGPLTT1Qyuu7rr7/O7XBFAWVRgvTVV19x+vRpLl26xPTp0ylcuDAzZ86kcuXKlNTPjCqsTGqQhBBPEW/9UiOgS5A8PT0pUqSIeqxTp07UrFmTQ0BjIF5/LuDhQ3YChql69+7dm1chiwLmiWqQfH19KVKkCL6+vvj4+ODk5JTh2jzCCmSYvxDiKeLi4qJuFylSBI1Gw8KFC+nQoQM7duzAycmJmjVrAnAOmNimDYp+3x/4C2gLXLx4MY8jFwWFRQnSu+++S6NGjShSpAhjx44lKSmJsWPHcu3aNSIjI60dowAZ5i+EeKpMmDBB3f7kk08AeOGFF1i7di1NmjQBICwsTJ2J+5/Tp7n+009s199TGPgTeObOHRIS1EoloqOjGTx4MBs3bsyDdyGeZhbNg/TJJ59QtGhRJk6cyPPPP0/FihWtHZcwJ6PYhBBPkdatW/Prr7/i4OBAq1atMrymUKFCVK5cmePHjxMZGUlow4Y8AH4Angdc9NvXPvkEL/1C6e+88w5r1qzhiy++IC4uLt2CukJkl0UtSJGRkYwfP559+/bRuHFjSpQowcsvv8yiRYtkxtPcIqPYhBBPEY1GQ7du3ejatWuW1xm62QDu3LnDI3RTAHyjXyjdAQiaOhUmTgRFMZmMcvny5VaPWxQcFiVINWrUYNiwYfz666/ExcWxbt06XFxcGDx4MFWqVLF2jAJkFJsQokAyTpAMtMCa9u352PjgpElo33rL5ENtyZIlMk+SsJhFCZKiKERERDBr1iy6dOlCy5Yt+e6776hWrRrDhg2zdowCZBSbEKJA6tixY4bLibRq3Zr3gBFGxxwWLuQHdF1vAMeOHSMqKir3gxRPJYsSJD8/P+rXr88PP/xAhQoVWLZsGTdv3iQiIoLZs2dbO0YBaoKkBR4iLUhCiIIhNDSUI0eOcPDgQZYtW4aLiwtNmjRRJ5acC8yrVw9FP/FkD3TF2x76+y9cuGCDqMXTwKIi7e+++46mTZtKK0Ze0tcgGdqR5HsvhCgoDOu01a5dm27dulG4cGFSUlJwcHBAq9Wy4N49Ivz8+OLGDdzRDf//C+gAxMXFAZCYmIiHh0dmLyFEOha1IHXs2BGtVsvMmTMZMGAAAwYMYNasWcTHx1s7PmGgb0GSBEkIUZB5enqi0WhwcXEhKCgIgJMnT7Lsxg3aAnf014UDO4BHZ87Qr18/vLy80i1uu2XLFpkOQGRKo1hQwXbgwAHatWtHoUKFCA8PB2D//v08fPiQTZs2Ubt2basHam8SEhLw9vYmPj4+b5KVokXh5k3OAiFAWlqaOj+IEEIURE2bNmXnzp0mx8KAjUCQfv9O4cI0un+fk/p9w0fe3r17adCgAQCrV6/mueeey4uQhR3I7ue3RZ+wI0eOpEuXLly4cIFff/2VX3/9lfPnz9OpU6cMV2kWVmDUxebp6SnJkRCiwGvcuHG6Y1FAI8CwMqjv/fvs4L+lSbRaLQC///67ek+3bt1yM0yRT1n0KXvgwAHGjBljshqzk5MT//vf/zhw4IDVghN6aWmQlATIEH8hhDDo2LFjhscvolu/LUK/b1iapA1w+fJlgHTD/0+dOpVLUYr8yqIEycvLi5iYmHTHL126JB/eueHBA3VThvgLIYROw4YNMz1+U6OhBbrECHRLk6wF4hctAuD69esm96xatSq3whT5lEUJUo8ePejfvz8rVqzg0qVLXLp0iZ9++okBAwbQs2dPa8coZKFaIYRIx8nJiXLlyqn7VatWpXv37sydOxc/Pz/uoRvJZkh9XICwKVPgiy/SJUgZ/dEvCjaLhvnPmDEDjUZDnz59SE1NBcDZ2ZlBgwapiw4KK5KFaoUQIkNLliyhZcuWuLm58eeff1KmTBkAihYtyq1bt3gEvAgsBAYCGkWBwYPpUrw4a42eExsbm+exC/tmUYLk4uLC3LlzmTp1KmfPngWgfPnyuLu7WzU4oScL1QohRIZatGjBmTNn8PDwoHjx4urxokWLcvKkbuyaFngdiAPe1Z9/PTaWZGAYoCAJkkjPogTJwN3dnbCwMEC38KDIJbLMiBBCZCokJCTdsaJFi6Y7Nh5dkmRY72EIugLuPkiCJNKzeKz4kiVLCAsLw83NDTc3N8LCwvjqq6+sGZtdmj9/PqGhodSrV+/xF1uLUQ2SjGITQojHyyhBApiDLiFK1e+/BPwB3IuNVacAEAIsTJAmTJjA8OHD6dy5MytXrmTlypV07tyZkSNHMmHCBGvHaFcGDx7M8ePH2b9/f969qLQgCSFEjvj7+2d6bjnwHLp1LQHaARvS0rh1+jS7du3ivffe48qVK7kfpLBrFnWxLViwgMWLF5uMWOvSpQvVq1dn6NChTJo0yWoBCtIlSCWkBUkIIbKUUQvSxx9/zPjx4wHdkP+26FqPfIEGwMMOHXjr7l2O3LnDnj172LJlSx5GLOyNRS1IKSkp1K1bN93xOnXqqKPahBWZDfOXFiQhhMhaYGCgyb6npyclSpQwObYLeLZQIa7q9wudP88fd+5QCdi6dSt79uyRz7QCzKIEqXfv3ixYsCDd8UWLFvHKK688cVDCjAzzF0KIHDGss2ZQpEiRDFuVXOrUoTH/LU1SCtiJbmmShg0b8r///S+XIxX26omLtAcMGMCAAQOoVq0aixcvxsHBgVGjRqlfwgpkmL8QQuRI6dKlKV++vLrv5uZGQEBAuuvq16/PBaAJcMrDAzBdmmTevHmkpaXx6NEjTpw4kW6JEvH0sihBioqKonbt2hQtWpSzZ89y9uxZ/P39qV27NlFRUURGRhIZGcmhQ4esHG4BJV1sQgiRY8ajjU+ePJlhC1KTJk0AuAG00GrZrj9uWJrk+bQ0jh49Svv27QkNDWXcuHG5HbawExYVaW/bts3acYismHWxSQuSEEI8XqNGjfjpp58A8PPzS5cgeXh4ULVqVXX/2sOHPAv8CHRDtzTJT8DyN95g2759AEybNi3DFSOSk5NZuXIl5cqVy3SNuMwoisKBAwcIDg5OVzslbMfiLjaRh2SYvxBC5Fi/fv0IDg4GYPbs2bi7u+Pq6qqeL1++PEWKFDG55xHQHTjRqBGg+5Dsu28fkx/zWl9++SW9evWiUaNGXL169TFXm/rpp58IDw+ncuXKJCQk5OhekXskQcoPpAZJCCFyzNPTkyNHjnD8+HH69OkDwKNHj9Tz5cuXx9vbO919aYDLsmXMdnNTj70HfIWu2yUpKSndPcOGDVO3f/755xzF+fLLLwMQHx/PDz/8kKN7Re6RBCk/MJtJW1qQhBAie3x8fKhSpUqG58qWLYujoyM+Pj7pzhUPCuLvdu0Ygm4tN4D+wBqgV7du1KtXj336bjdzTk6Wr+KVkpJi8b3CuiRByg+MWpAeoOs3F0II8WR8fX0BXX2SMU9PT9zd3WnevDnzgRfRdb0BdATe2bCBCwcOMGXKlAyf6+zsnGsxi7yT4wQpJSWF1q1bc+bMmcdfLKxDnyClAG5eXrIwsBBCWIFh2L95gmQolG7RogUAq9AtR3JXf74+ukkmb2ay5JSjo6O6ff/+ffr168eIESNkrbd8JscJkrOzM0eOHMmNWERm9F1sUn8khBBP5rvvvgN0E0e+9NJLQOYJUvXq1dXtv4GmgGGFtorAL1evkrxnT7rXePjwobo9ceJEli1bxty5c1mxYsVj45N5luyHRV1svXr1YsmSJdaORWRG34Ik9UdCCPFkXn75ZSIiIjh+/Lj6+9Q8QQoKCgJ0LUG//PILw4YNY86cOUQBDYHj+usCAcfWrWHzZpP7E43KIr744gt1Oztru0mCZD8sqiRLTU1l6dKlbNmyhTp16qSriZk1a5ZVghM6SmIiGqQFSQghnpRGo6FWrVomxwy1SAbGM3A3btyYxo0bqz0nl9DNuv27/r+ODx6gbd+elwHD+LP7RgNrjEe8Zad+1HiUnbAtixIkw0zaAKdPnzY5J/UxucCoi01akIQQwrrMW5DKlSuX7prixYur23eAtugSom6AQ1oa3wPFgZmYtiAZc3d3f2wsDx48yG7YIpfJTNr2LjkZjX41aeliE0II68tOglSkSBGCg4O5dOkSVatW5dixY7wAfAa8pb9mBlACOGnUgmTMuDYpM5Ig2Q8Z5m/vZJJIIYTIVdlJkBwcHFi3bh2zZs3ir7/+4tVXX0ULDAbGG103Eui3aRM8epRu1Nrdu3cfG4skSPbD8tmsgOPHjxMTE0NycrLJ8S5dujxRUMKILDMihBC5yjxBKlmyZIbXhYWFERYWBsDSpUspVaoUH374IVOAa8CX6D5UG8bEcK12bV41m6VbEqT8xaIE6dy5c3Tr1o2jR4+i0WjUqntD/VFaWpr1IizozGbRlhYkIYSwLvMibeN5jLLSr18/pk2bRlJSEkvRJUkrAXcg8PhxPgEOA7H666Oioli6dCldunTB398fSD9qTRIk+2FRF9vw4cMpW7YsN27cwN3dnWPHjvHPP/9Qt25dtm/fbuUQCzhpQRJCiFxl6dIgZcqU4cSJE+r+OqAlcFO/XwP4F6ik3z937hz9+/enX79+6j3m67pJgmQ/LEqQdu/ezaRJk/D398fBwQEHBweaNGnC1KlTTRbsE1YgNUhCCJGrateurU4IOXPmzBzdW6ZMGZP9fUAj4LzhPLpZtxsYXbN27Vpu3tSlUeYJkSRI9sOiBCktLU39oPb39+fq1asAlC5dmlOnTlkvOmHSxSYtSEIIYX2urq5ERkaybds2hg8fnuP73377bZP9M+gmlIzQ7xcBtgIdjK5Zt24dCxYsMJk+ALI30k3kDYvaFcPCwjh8+DBly5alfv36TJ8+HRcXFxYtWpRh9b94AkYtSFKDJIQQuSMwMFBtRcqpSZMmkZCQwOLFi9Vj14EW6NZxa4uuLuk3oC+6+ZP++usvli1blu5Z0oJkPyxqQXrvvffU4YuTJk3i/PnzNG3alHXr1jFv3jyrBljgSQ2SEELYNXd3dxYtWpTu+D2g6vnz/K1PvJyA74EhwLVr1zJ8liRI9sOiFqR27dqp2yEhIZw8eZLbt2/j6+srM2lbm1kXm7QgCSFE/lG8dGlGN2vG8Z9/ZpD+2GfABxs3sjGD6yVBsh9WmyjSz89PkqPcYNbFJi1IQgiRP/To0QONRkOJUqV4C5hsdO4DYB5g/qmZnbmSRN7IdgvSqFGjsv1QWazWimQUmxBC5Dvvv/8+77//PgDNmjVjxowZTABuAXP01wxFV8DdF0jVH7t//z5r1qzhueeey9uARTrZTpAiIyNzMw6RGRnFJoQQ+UKbNm3YsmULAAMHDsTZ2RmApk2bqtfMBW4DS9F9AL8M+AAvAIbxa5MmTaJDhw64uLjkVegiA9lOkGSBWhsxakF6oNFkazVoIYQQeW/x4sVMnDiRFi1aEBwcrB738fHBxcVFXZZrOXAX+BlwQzf8fxPQCYhH1yAREhLCsWPHpNfAhiwq0p40aVKm5zQajdqsKKzAKEHCw0PqvIQQwk6VKVMmw6H7ADt27GDSpEls3LiR1NRU/gDaAb8D3kAT4F8nJ1qnpnINuHTpEgsWLOB///tfnsUvTGkU84VgsqFWrVom+ykpKZw/fx4nJyfKly9PREREJnc+PRISEvD29iY+Pj53u726dYM1awCoGxTEgStXcu+1hBBC5KohQ4Ywf/58db8msBEI0O+fRTdv0nkgPDyc69evU7VqVf744w8cHKw2rqpAy+7nt0UtSBnVIyUkJNCvXz+6detmySNFZoxqkDTS1CqEEPlaamqqyf4hIPqbbwiYMAFiYiiPbmmSZ4B9+/YBcPHiRXbt2mVSyyRyn9XSUS8vLz788EPpXrMyxShBcvHxsV0gQgghntijR4/SHXOrVg127eKSvjWjOPAPuuVKDJo1a8bQoUOpW7cumzZtypNYCzqrttfFx8cTHx9vzUcWeNp79wDd6AYPGcEmhBD5WkbLcXl5eUHJkszp1o29+mO+wBZ0dUoGn3/+OQcPHuS9994DIC4ujmXLlhEbG5vbYRdIFnWxmS8noigKsbGxLF++nPbt21slMKGj6BMkGeIvhBD534gRI/j+++9NFnYvVqwYAM6BgbQGVvPf+m1/AL2BFUbP2L9/PwCvvPIKmzdvpm7duuzbt08G8ViZRQnS7NmzTfYdHBwoWrQoffv2Zdy4cVYJTOjpR7HJQrVCCJH/eXp6cvz4caKjo5k+fTrt27dXf7drNBoS0Q33/w7oDjijW9zWF1ho9Bxvb28SEhIAOHDgAFFRUVSrVi0v38pTz6IE6fz589aOQ2RC81A3dZi0IAkhxNPBwcGBihUr8tVXX5kcNxRwJwMvAYWDg2l/6RIOwALAH/hIf60hOTL4+eefJUGyMhkzaM8UBQejBElakIQQ4ulVtGhRdbtMuXK0v3iRqUbnJwOzSb9+G8DKlSuxYNYekQWLWpAyW5dNo9Hg5uZGSEgIXbt2xc/P74mCK/AePkSj/wcvC9UKIcTT7c0332T27Nncvn2bH374ATQaXGfOZPTbbzNDf80IwA/oz3/rtwGcOnUKBwcHfHx8GDduHEOHDqVQoUJ5/RaeKhZNFNmyZUsiIiJIS0ujUqVKAJw+fRpHR0cqV67MqVOn0Gg07Ny5k9DQUKsHbS0zZszg66+/RqPRMHbsWHr16pXte/Nkosi4OAjQTR/2B3D5iy8YNGhQ7ryWEEIIm0tKSuLhw4f4+vqqx+rUqUP1iAi+Ahz1x34HegBJmTzn888/Z/DgwWqrkhRw/ye7n98WdbF17dqVNm3acPXqVQ4ePMjBgwe5fPkybdu2pWfPnly5coVmzZoxcuRIi99Abjt69Cg//PADBw8eZP/+/Xz++efcvXvX1mGZMlpmRGqQhBDi6efm5maSHBl8g25BW8MsSl3QzcDtBRle/+OPP3L27FmCg4OpUaMGicbLVolssShB+vTTT5k8ebLJB7a3tzcffPAB06dPx93dnQkTJnDw4EGrBWptJ06coGHDhri5uVGoUCFq1KjBhg0bbB2WKaNJImUUmxBCFEzOzs4ArAGeBe7pjzcD/gZqlyiR7p6AgABeeuklrly5wtGjR/nxxx/zJtiniEUJUnx8PDdu3Eh3PC4uTq2s9/HxUVcutsQ///xD586dCQoKQqPRsEa/Hpmx+fPnU6ZMGdzc3Khfv746LXt2hIWFsX37du7evcudO3fYvn07V+xtnTNpQRJCiAJvxowZ6va7mzbxde/exOn3awJfR0dT2uye+Ph4Dhw4oO5n9JktsmZRkXbXrl157bXXmDlzJvXq1QN0E1eNHj2a5557DtCtIVOxYkWLA0tMTKRGjRq89tprPP/88+nOr1ixglGjRrFw4ULq16/PnDlzaNeuHadOnSJAX7dTs2bNdOveAGzatInQ0FCGDRtGq1at8Pb2pkGDBjg6Oqa71qYkQRJCiAKvSZMm/P333zg7O9OwYUOioqJounw5m4BSQHBSEn8DLdEtcgtw9uxZk2ekpKTkbdBPA8UC9+7dUwYMGKC4uLgoDg4OioODg+Li4qIMHDhQuX//vqIoihIZGalERkZa8vh0AGX16tUmx8LDw5XBgwer+2lpaUpQUJAydepUi16jf//+yp9//pnp+aSkJCU+Pl79unTpkgIo8fHxFr1etvz2m6KAooDyLiinT5/OvdcSQgiRL+zYsUMBlJKgXC5cWP2ciAGlLChk8DVo0CBbh2034uPjs/X5bVEXW+HChVm8eDG3bt0iMjKSyMhIbt26xaJFi/Dw8AB0rTc1a9a0QgqXXnJyMgcPHqRNmzbqMQcHB9q0acPu3buz/RxDk+OpU6fYt28f7dq1y/TaqVOn4u3trX4FBwdb/gayy6gFSYb5CyGEAF2L0oQJE2jwwgs47dzJFf1nQzCwHt00AOakiy3nLOpiMyhcuDDVq1e3VizZdvPmTdLS0tT1awyKFSvGyZMns/2crl27Eh8fj4eHB19//TVOTpl/O8aNG2cy/1NCQkLuJ0lmXWxSpC2EEALgww8/VLdHde7MgO+/JxSohG4KgDaYTgGwatUqYmJiKFWqVN4Gmo9ZnCBt3bqVrVu3cuPGDbRarcm5pUuXPnFgeSEnrU2urq64urrmYjQZMBrF9lCjkUm/hBBCpJNWpAjtgT1AcaAxsBzdPEnGn85LliwxSaxE1izqYvvwww955pln2Lp1Kzdv3uTOnTsmX7nN398fR0dHrl+/bnL8+vXrBAYG5vrr5xmjFiStu7tM9CWEECIdd3d3YoCO6MoxQDdn0gyz6w4fPpynceV3FrUgLVy4kG+++YbevXtbO55scXFxoU6dOmzdulUdNafVatm6dStDhgyxSUy5wnhiL31tlxBCCGHM3d0dgEh0idFajQZHRWEkcBGYq7/OxcXFNgHmUxa1ICUnJ9OoUSNrx2Li/v37HDp0iEOHDgFw/vx5Dh06RExMDKBbD27x4sUsW7aMEydOMGjQIBITE3n11VdzNa48JQmSEEKIxzAkSKCbXXuX0bJZc4BX9NuG1SIURSEuLg6RNYsSpAEDBugW0stFBw4coFatWtSqVQvQJUS1atViwoQJAPTo0YMZM2YwYcIEatasyaFDh9iwYUO6wu38THvvnrrtICPYhBBCZMC8ZajImDFMMtr/Bt3SJIYEqWPHjhQrVoxBgwbRsWNHvvrqK5P7Y2JiMpxDsKCxqIstKSmJRYsWsWXLFqpXr65Og24wa9asJw6sRYsW6iJ7mRkyZMjT1aVmJvXuXQz/7B0lQRJCCJENVatWpYaDA0W0Wgaj+6D/GXjt0iU2bdrE+vXrAV25DMC6det4+eWXcXd358svv+TNN9+kfv367N69u0DXvlqUIB05ckSd4ygqKsrkXEH+Zlpbmn7ZFgBnHx/bBSKEECJfOXDwIGPeeYc28fFU2r8fV+DLa9doncl8fwkJCbi7u/Pmm28CsHfvXqKjo6lQoUIeRm1fLEqQtm3bZu04RAa0kiAJIYSwQM2aNdm4eTOkpLC9aFFaxMdTGN1Eks2AY2bX3zeaVsbg0aNHeRCp/XqiiSKPHz9OTEyMyaK0Go2Gzp07P3FgAhT9P1gt4Obra9tghBBC5D/OzkytUQPtP//QCt0s21uAVsAJo8skQUrPogTp3LlzdOvWjaNHj6LRaNRaIUP3WlpamvUiLMj0o9gSAS9vb9vGIoQQwi7VqFFD3W7ZsmW68x5FitAVXWJUHwgE/sI0STpz5gwODqbjtjJKmgoSi0axDR8+nLJly3Ljxg3c3d05duwY//zzD3Xr1mX79u1WDrHg0jx4AMgyI0IIITLXrFkzxo0bR6dOnfjmm2/SnS9atCj3gWeBA/pjhiSpin7/xRdfNEm0AO4ZjaQuiCxqQdq9ezd//fUX/v7+ODg44ODgQJMmTZg6dSrDhg0jMjLS2nEWSI5JupV0ZKFaIYQQWZkyZUqm58LCwgC4C7QFNgN1ybglydjly5d57bXXCAgIoH79+syYMYMRI0bQvXt3K0dvnyxKkNLS0tQWDX9/f65evUqlSpUoXbo0p06dsmqABZmTvv9XWpCEEEJYyjDqHDJOkrahS5KOm903fPhwkxpjgH///fexU/A8LSzqYgsLC1PXdKlfvz7Tp09n165dTJo0iXLlylk1wAIrLQ2nlBRAX4MkLUhCCCEsUL16dZP9u+iSJEN3WzHgH3QJkzHz5Mjg1q1bVo3PXlmUIL333ntotbo1gidNmsT58+dp2rQp69atY968eVYNsMAyWmbkPtKCJIQQwjLe3t5Mnz7d5Nhd4Nc332Sffr8Iuu62Ftl43okTGXXIPX0sSpDatWvH888/D0BISAgnT57k5s2b3Lhxg1atWlk1wALLKEGSFiQhhBBP4p133uH4cdNONJ+yZWmDrosNwBPdPEldHvOskydPWj9AO2RRgpQRPz8/mUXbmswSJGlBEkII8SSqVKlisl+iRAnuAR2A3/XH3IBVQC8yZ1g0/mlntQRJWJnR/BMyik0IIYQ1HDp0iDp16jBs2DDKly8PQBLwf8B3+mucgOXAyEyeMXny5Kd6HVQDjVJQytGtLCEhAW9vb+Lj43Mnedm1C5o0AWAm8NaDBxQqVMj6ryOEEKJAUhSFli1b8vfffwOgAeYBxqnPPHSJkjaD+69evUrx4sVzPU5ry+7nt7Qg2SujLraHGg1ubm42DEYIIcTTRqPRsG3bNq5du8ZLL71EEX9/hgITja4Zhq7LLaM/z9euXZsncdqKJEj2yqiLLdXNTeq7hBBCWJ1Go6FYsWL8+OOPXL9+nZ49ezIJ6Aek6K95DtgOBJjdu2/fPp5mkiDZK6MWpDTpWhNCCJHLHBwc+OGHHwgJCWEZ0B6I158LB3YDoUbXL168mKNHj+Z1mHlGEiR7ZZQgKe7uNgxECCFEQfLrr79SvXp1tgJNgEv64+WAPZhOAzB58uS8Di/PSIJkp9KMFwmUBEkIIUQeqVatGocPH2b79u1c9PSkARChP+cJ/Aa8q9+/fPkyt2/fVu+Njo7mueeeY+rUqfl+SRJJkOxUstFU7hqZA0kIIUQea968OQkJCVxRFGonJkKPHuq5j4GfgMO7d1OkSBF1JNzs2bP57bffePfdd/ntt99sE7iVSIJkp5Lv3lW3HSRBEkIIYUvu7vDjjzBlCop+0FAPYCdQBmjRogUffPABX3zxhXrL0qVLbRGp1UiCZKdSjRIkJ29v2wUihBBCAGg0MG4cyurVJOgP1ULX/dYZ+PDDD00u/+OPP/joo4/YvXs3ixYt4p5x6Ug+IAmSnUqLj1e3HSVBEkIIYSccunalAXBGv++LbqmSTwBHs2vff/99GjVqxBtvvMHEiRPJTyRBslNao3mQXP38bBiJEEIIYapshw7UBX4xOjYG+AvIbG7t2bNn53pc1iQJkp1SjBIkF19fG0YihBBCmFq6dCkJQHdgOP9NKtkMiATa2iowK5IEyU5pjOZBKuTvb8NIhBBCCFPFihXj/v37rF+/np67d9OM/+ZLKgZsQreOqIvZfb/88gvm4uPjmTx5Mn/88UeuxpxTTrYOQGTM4cEDAJIBD2lBEkIIYWc8PDx49tlnuX79OnvQFWwvRzcDN8AooDXwMnBcf6x79+7Mnj2b2NhYatasycmTJzl16hQrVqwA4OLFi5QqVSqP30nGJEGyU45JSQAkQparDQshhBC25Kevk70FdASGAtMBV6AGcAAYDRgmABg5cmSmz9q+fTt9+vTJxWizT7rY7JTjo0eALkHylHmQhBBC2ClnZ2datmwJQPGgIE4/+ywjGjXimP58IWA+8CcQ9JhnJScn516gOSQJkp1y0f8juY+0IAkhhLBv69atY8+ePcTExLB+/XoW7NrFqrFjmWd0TUfgGNAvi+dIgiQeyyVFNyZAWpCEEELYOzc3N+rXr4+j438zIRUuWpThQAcgVn/MB/gaWAeUzOA5t2/fZsSIEbRr185kjTdbkBoke5ScjJN+kb9EoJS0IAkhhMhnwsLCAFgPVAVmA33159oDUcDbwBKje95//311u0iRIiQlJeHq6poX4aYjLUj2yGgOJOliE0IIkR81atRI3XYqWpR+6LrZruiPeQNfoZtcslImz1i/fn0uRpg1SZDskdEcSA81Gptlz0IIIYSlChcuzKJFi2jUqBE//fQTLVu2ZB261iTjZWxbAkeASYCb2TMePnyYR9GmJwmSPTJKkJJdzKfZEkIIIfKHgQMHsmvXLlq1asWqVav4+eef2X/6NP2BZ4Fz+utcgPfRdbs9Y3R/WlpaHkf8H0mQ7JFRF1uqJEhCCCGeAr6+vnTv3l2dCHIjEAZ8X6YMhrFr5fXHVwAlgEf6KW9sQRIke2TUgpTqZt7gKIQQQuRfxmUjD4FXzp/n+saNbDe65kXgFFBp1SrQavM2QD1JkOxQany8up0mCZIQQoinXMm2bRlerRp9gBv6Yx5AkbNnwcE2qYokSHYo6eZNdVvx8LBhJEIIIYT1TZkyhcKFCzN37lwANBoNm7dsoc2yZVQC5qEbxb2ifn2bxSgJkh1KunXrvx1JkIQQQjxlxo0bR3x8PMOGDVOPBQQE8Morr3AXGA6UBqp07GijCCVBskvJd+6o25rChW0YiRBCCJE7HDLoOnN0dGTu3LnUrl2bFZs306NHDxtEpiMzadsh4wTJUSaJFEIIUYAMGzbMpGXJVqQFyQ4ZF2k7envbMBIhhBCiYJIEyQ6lJSSo284+PrYLRAghhCigJEGyQ8q9e+q2i6+vDSMRQgghCiZJkOyQYjSTtmuRIjaMRAghhCiYJEGyQxqjmbTdJEESQggh8pwkSHbIwWj1Yo+AABtGIoQQQhRMkiDZIcekJHXb3d/fhpEIIYQQBZMkSHbIWb968QPAS4q0hRBCiDwnCZIdck5JASAR8PT0tG0wQgghRAEkCZIdcpEESQghhLApSZDsUKG0NAASNRpcXV1tHI0QQghR8EiCZG8UBTetFoAkR0cbByOEEEIUTJIg2ZuHD9X/KY+cZC1hIYQQwhYkQbI3RrNoJ7u42DAQIYQQouCSBMnOJN+5o26nSoIkhBBC2IQkSHbmQVycup3m5mbDSIQQQoiCSxIkO/Pgxg11W+vubsNIhBBCiIJLEiQ7k3TrlrqtSIIkhBBC2IQkSHbGOEGicGHbBSKEEEIUYJIg2ZlHt2+r2w6SIAkhhBA2IQmSnUm9e1fddvDysl0gQgghRAEmCZKdSY2PV7edvL1tGIkQQghRcEmCZGfS7t1Tt519fW0YiRBCCFFwSYJkZ5SEBHXbRRIkIYQQwiYkQbI3iYnqpluRIjYMRAghhCi4JEGyM5oHD9RtSZCEEEII25AEyc44PnyobnsEBNgwEiGEEKLgkgTJzjgmJanb7pIgCSGEEDYhCZKdcUpOVrcLFytmw0iEEEKIgksSJDvjok+QtICLTBQphBBC2IQkSHbGNTUVgESNBjQaG0cjhBBCFEySINkZt7Q0AB46yP8aIYQQwlbkU9jOFNJqAUhydLRxJEIIIUTBVWATpPnz51OmTBnc3NyoX78++/bts3VIPHr0iMKGbWdnm8YihBBCFGQFMkFasWIFo0aNYuLEiURERFCjRg3atWvHjRs3bBpXwp07uOm3U1xcbBqLEEIIUZAVyARp1qxZDBw4kFdffZXQ0FAWLlyIu7s7S5cutWlc969fV7dTXV1tGIkQQghRsBW4BCk5OZmDBw/Spk0b9ZiDgwNt2rRh9+7dmd736NEjEhISTL6sLfXuXXVb6+aW+YVCCCGEyFUFLkG6efMmaWlpFDObhLFYsWJcu3Yt0/umTp2Kt7e3+hUcHGz12CoEBanbNZo0sfrzhRBCCJE9BS5BstS4ceOIj49Xvy5dumT9F3F1hW7d4Jln0NSoYf3nCyGEECJbnGwdQF7z9/fH0dGR60b1PgDXr18nMDAw0/tcXV1xze26oFKl4Ndfc/c1hBBCCPFYBa4FycXFhTp16rB161b1mFarZevWrTRs2NCGkQkhhBDCXhS4FiSAUaNG0bdvX+rWrUt4eDhz5swhMTGRV1991dahCSGEEMIOFMgEqUePHsTFxTFhwgSuXbtGzZo12bBhQ7rCbSGEEEIUTBpFURRbB5EfJSQk4O3tTXx8PF5eXrYORwghhBDZkN3P7wJXgySEEEII8TiSIAkhhBBCmJEESQghhBDCjCRIQgghhBBmJEESQgghhDAjCZIQQgghhBlJkIQQQgghzEiCJIQQQghhRhIkIYQQQggzBXKpEWswTECekJBg40iEEEIIkV2Gz+3HLSQiCZKF7t27B0BwcLCNIxFCCCFETt27dw9vb+9Mz8tabBbSarVcvXoVT09PNBqN1Z6bkJBAcHAwly5dkjXeCiD5/5+xevXqsX//fluHkSlbxJdbr2mt5z7pcyy5P6f3ZOd6+ZnMXH79uVQUhXv37hEUFISDQ+aVRtKCZCEHBwdKliyZa8/38vKSH8YCTP7/m3J0dLTr74ct4sut17TWc5/0OZbcn9N7cnK9/Eyml59/LrNqOTKQIm0hhN0bPHiwrUPIki3iy63XtNZzn/Q5ltyf03vs/d+VvbP379+TxiddbHYmISEBb29v4uPj7TozF7lD/v8LYV/kZ7LgkhYkO+Pq6srEiRNxdXW1dSjCBuT/vxD2RX4mCy5pQRJCCCGEMCMtSEIIIYQQZiRBEkIIIYQwIwmSEEIIIYQZSZCEEEIIIcxIgmRH5s+fT5kyZXBzc6N+/frs27fP1iEJIYQQBZIkSHZixYoVjBo1iokTJxIREUGNGjVo164dN27csHVowk5069YNX19fXnjhBVuHIoQALl26RIsWLQgNDaV69eqsXLnS1iEJK5Jh/naifv361KtXj88//xzQrfUWHBzM0KFDGTt2rI2jE/Zg+/bt3Lt3j2XLlvHLL7/YOhwhCrzY2FiuX79OzZo1uXbtGnXq1OH06dN4eHjYOjRhBdKCZAeSk5M5ePAgbdq0UY85ODjQpk0bdu/ebcPIhD1p0aIFnp6etg5DCKFXvHhxatasCUBgYCD+/v7cvn3btkEJq5EEyQ7cvHmTtLQ0ihUrZnK8WLFiXLt2zUZRCWv6559/6Ny5M0FBQWg0GtasWZPuGqlBEyJvWfPn8uDBg6SlpREcHJzLUYu8IgmSEHkgMTGRGjVqMH/+/AzPSw2aEHnPWj+Xt2/fpk+fPixatCgvwhZ5RBIkO+Dv74+joyPXr183OX79+nUCAwNtFJWwpvbt2/PRRx/RrVu3DM/PmjWLgQMH8uqrrxIaGsrChQtxd3dn6dKleRypEAWHNX4uHz16xHPPPcfYsWNp1KhRXoUu8oAkSHbAxcWFOnXqsHXrVvWYVqtl69atNGzY0IaRibwgNWhC2J/s/FwqikK/fv1o1aoVvXv3tlWoIpdIgmQnRo0axeLFi1m2bBknTpxg0KBBJCYm8uqrr9o6NJHLsluD1qZNG7p37866desoWbKkJE9C5KLs/Fzu2rWLFStWsGbNGmrWrEnNmjU5evSoLcIVucDJ1gEInR49ehAXF8eECRO4du0aNWvWZMOGDel+OEXBtWXLFluHIIQw0qRJE7Rara3DELlEEiQ7MmTIEIYMGWLrMEQekxo0IeyP/FwK6WITwsakBk0I+yM/l0JakITIA/fv3yc6OlrdP3/+PIcOHcLPz49SpUoxatQo+vbtS926dQkPD2fOnDlSgyZELpOfS5EVWWpEiDywfft2WrZsme543759+eabbwD4/PPP+fTTT9UatHnz5lG/fv08jlSIgkN+LkVWJEESQgghhDAjNUhCCCGEEGYkQRJCCCGEMCMJkhBCCCGEGUmQhBBCCCHMSIIkhBBCCGFGEiQhhBBCCDOSIAkhhBBCmJEESQghhBDCjCRIQgghhBBmJEESQgghhDAjCZIQIte0aNGCESNG2DqMJ5bf30dexz969Giee+65PHs9IXKDJEhCCCGs6tChQ9SsWdPWYQjxRCRBEkLYteTkZFuHIDKQ1f+Xw4cPS4Ik8j1JkIR4Cm3YsIEmTZrg4+NDkSJF6NSpE2fPnlXPt2jRgmHDhvG///0PPz8/AgMD+eCDD0yece/ePV555RU8PDwoXrw4s2fPTtdVU6ZMGebMmWNyX82aNdM9K7txGWIbMmQII0aMwN/fn3bt2mX4rBYtWjB06FBGjBiBr68vxYoVY/HixSQmJvLqq6/i6elJSEgI69evV+959OgRw4YNIyAgADc3N5o0acL+/ftNnpuYmEifPn0oXLgwxYsXZ+bMmeleW6vVMnXqVMqWLUuhQoWoUaMGv/zyS4ZxGseb1fc8O99LS96zQWpqKkOGDMHb2xt/f3/ef/99FEXJ9nvK7v+Xy5cvc/PmTWrUqJHl90MIeycJkhBPocTEREaNGsWBAwfYunUrDg4OdOvWDa1Wq16zbNkyPDw82Lt3L9OnT2fSpEls3rxZPT9q1Ch27drF77//zubNm9mxYwcRERG5HpchNhcXF3bt2sXChQszfd6yZcvw9/dn3759DB06lEGDBtG9e3caNWpEREQEzzzzDL179+bBgwcA/O9//2PVqlUsW7aMiIgIQkJCaNeuHbdv31af+c477/D333/z22+/sWnTJrZv357ufU+dOpVvv/2WhQsXcuzYMUaOHEmvXr34+++/s3z/j/ueZ0dO37PxfU5OTuzbt4+5c+cya9Ysvvrqqxy9p+z8fzl06BDe3t6ULVs2R+9LCLujCCGeenFxcQqgHD16VFEURWnevLnSpEkTk2vq1aunjBkzRlEURUlISFCcnZ2VlStXqufv3r2ruLu7K8OHD1ePlS5dWpk9e7bJc2rUqKFMnDhRfR3j6x8Xl+GeWrVqPfY9mb+H1NRUxcPDQ+ndu7d6LDY2VgGU3bt3K/fv31ecnZ2V77//Xj2fnJysBAUFKdOnT1cURVHu3bunuLi4KD///LN6za1bt5RChQqp7yMpKUlxd3dX/v33X5N4+vfvr/Ts2TPb8SqK6ff8cd9LS96z8X1VqlRRtFqtemzMmDFKlSpVsv2esvv/ZfLkyUqzZs0ee50Q9k5akIR4Cp05c4aePXtSrlw5vLy8KFOmDAAxMTHqNdWrVze5p3jx4ty4cQOAc+fOkZKSQnh4uHre29ubSpUq5XpcAHXq1MnW84zfg6OjI0WKFKFatWrqsWLFigFw48YNzp49S0pKCo0bN1bPOzs7Ex4ezokTJwA4e/YsycnJ1K9fX73Gz8/P5H1HR0fz4MED2rZtS+HChdWvb7/9Nl13YVbxgun3PLty8p6NNWjQAI1Go+43bNiQM2fOkJaWlu33lJ3/L4cOHZLuNfFUcLJ1AEII6+vcuTOlS5dm8eLFBAUFodVqCQsLMymsdXZ2NrlHo9Gk6+p6HAcHB5M6FoCUlJQnigvAw8MjW6+f0XswPmZICHL6vrJy//59ANauXUuJEiVMzrm6umZ5b1bf8+x+L3PjPWf3PWXn/8uhQ4fo0KFDtl9bCHslLUhCPGVu3brFqVOneO+992jdujVVqlThzp07OXpGuXLlcHZ2Nilgjo+P5/Tp0ybXFS1alNjYWHU/ISGB8+fP51pcT6J8+fJq/YxBSkoK+/fvJzQ0VL3G2dmZvXv3qtfcuXPH5H2Hhobi6upKTEwMISEhJl/BwcEWx5eT76UljN8TwJ49e6hQoQKOjo5We0/37t3j3LlzMoJNPBWkBUmIp4yvry9FihRh0aJFFC9enJiYGMaOHZujZ3h6etK3b1/eeecd/Pz8CAgIYOLEiTg4OJh007Rq1YpvvvmGzp074+Pjw4QJE3B0dMy1uJ6Eh4cHgwYNUt9TqVKlmD59Og8ePKB///4AFC5cmP79+/POO+9QpEgRAgICGD9+PA4O//0t6enpyejRoxk5ciRarZYmTZoQHx/Prl278PLyom/fvhbFl5PvpSViYmIYNWoUb7zxBhEREXz22WfqCD1rvafDhw/j6OhI1apVrRa3ELYiCZIQTxkHBwd++uknhg0bRlhYGJUqVWLevHm0aNEiR8+ZNWsWb775Jp06dcLLy4v//e9/XLp0CTc3N/WacePGcf78eTp16oS3tzeTJ0/OtNXDWnE9iU8++QStVkvv3r25d+8edevWZePGjfj6+qrXfPrpp9y/f5/OnTvj6enJ22+/TXx8vMlzJk+eTNGiRZk6dSrnzp3Dx8eH2rVr8+6771ocW06+l5bo06cPDx8+JDw8HEdHR4YPH87rr7+unrfGezp06BCVK1d+bFejEPmBRjHv9BZCiAwkJiZSokQJZs6cqba4CCHE00pakIQQGYqMjOTkyZOEh4cTHx/PpEmTAOjatauNIxNCiNwnCZIQIlMzZszg1KlTuLi4UKdOHXbs2IG/v7+twxJCiFwnXWxCCCGEEGZkmL8QQgghhBlJkIQQQgghzEiCJIQQQghhRhIkIYQQQggzkiAJIYQQQpiRBEkIIYQQwowkSEIIIYQQZiRBEkIIIYQwIwmSEEIIIYQZSZCEEEIIIcxIgiSEEEIIYeb/AUCjROhF540EAAAAAElFTkSuQmCC", "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": ".venv", "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.12.7" } }, "nbformat": 4, "nbformat_minor": 0 }