bg01_friedmann_equations_distances.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Expansion of the universe and distances\n",
    "\n",
    "### Andrey Kravtsov\n",
    "\n",
    "#### The University of Chicago \n",
    "\n",
    "#### email: kravtsov@uchicago.edu\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Preliminaries: set up plotting parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# allow automatic reload of imported codes, if they are modified\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "# setup notebook for inline figures\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from codes.plot_utils import plot_pretty\n",
    "plot_pretty(dpi=225)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Expansion of the universe\n",
    "\n",
    "Solution of the Einstein's equations of general relativity in the Friedmann-Robertson-Walker (FRW) metric (homogeneous, isotropic space) gives the Friedmann equations for the expansion factor, defined as the distance between two freely expanding points and normalized to be unity at the current time. \n",
    "\n",
    "$$\n",
    "\\frac{\\ddot{a}}{a}=-\\frac{4\\pi G}{3}\\,(\\rho+3p) +\\frac{\\Lambda}{3},\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\left(\\frac{\\dot{a}}{a}\\right)\\equiv H(a)^2=\\frac{8\\pi G\\rho}{3} - \\frac{k}{a^2} + \\frac{\\Lambda}{3},\n",
    "$$\n",
    "\n",
    "where $H(a)$ is the Hubble expansion rate, equal to $H_0$ at $z=0$, $k$ is space curvature ($k = 0$ corresponds to a spatially flat, Euclidean universe, $k > 0$ to space with positive\n",
    "curvature, and $k < 0$ to space with negative curvature, such as saddle). $\\rho$ is the total energy density of the universe due to matter, radiation, and dark energy, and $p$ is the total pressure equal to the sum of pressures of each component. \n",
    "\n",
    "Here $\\Lambda$ is the cosmological constant that can be interpreted as the vacuum energy density $\\Lambda = 8\\pi G\\rho_{\\Lambda}=-8\\pi G p_{\\Lambda}$, i.e. equation of state of the vacuum energy is $p_{\\Lambda}=-\\rho_{\\Lambda}.$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The equations above are often recast with densities in units of the critical density of the universe \n",
    "\n",
    "$$\n",
    "\\rho_{\\rm crit}=\\frac{3H^2(a)}{8\\pi G};\\ \\ \\ \\Omega_{\\rm m}(a)=\\frac{\\bar{\\rho}_{\\rm m}(a)}{\\rho_{\\rm crit}(a)};\\ \\ \\ \\Omega_\\Lambda=\\frac{\\rho_\\Lambda}{\\rho_{\\rm crit}(a)};\\ \\ \\ \\Omega_k=1-\\Omega_{\\rm m}-\\Omega_\\Lambda.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "$$\n",
    "H(a)^2=\\frac{8\\pi G\\rho}{3} - \\frac{k}{a^2} + \\frac{\\Lambda}{3}=H_0^2\\left(\\Omega_{\\rm m0}\\,a^{-3} + \\Omega_k a^{-2} +\\Omega_{\\Lambda 0}\\right),\n",
    "$$\n",
    "\n",
    "A dimensionless Hubble function, $E(z)=H(z)/H_0$, is also often used:\n",
    "\n",
    "$$\n",
    "E(z)\\equiv H(z)/H_0=\\sqrt{\\Omega_{\\rm m,0}\\,(1+z)^3+\\Omega_k\\,(1+z)^2+\\Omega_{\\Lambda,0}}\n",
    "$$\n",
    "\n",
    "where redshift and expansion factor are related by $a=(1+z)^{-1}$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Thus, \n",
    "\n",
    "$$\\Omega_{\\rm m}(a)=\\frac{\\bar{\\rho}_{\\rm m}(a)}{\\rho_{\\rm crit}(a)}=\\frac{\\Omega_{\\rm m0}\\,a^{-3}}{\\Omega_{\\rm m0}\\,a^{-3} + \\Omega_k a^{-2} +\\Omega_{\\Lambda 0}}.$$\n",
    "\n",
    "The equation shows that at early epochs (decreasing $a$), $\\Omega_{\\rm m}(a)$ will approach $1$ regardless of the value of $\\Omega_{\\rm m0}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Distance - redshift relation\n",
    "\n",
    "When we observe an object with spectrum red-shifted by a cosmological redshift $z$, this redshift was accumulated during the time light travelled to us after it was emitted by the object. The emitted photons had to play catch up: in addition to covering the initial *physical* distance between us and the object  at the time light was emitted, the photons had to travel additional distance due to expansion of the universe between emission and observation. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The redshift of the photons during this travel was accumulated via a sequence of small redshifts $dz=H(z)/c\\delta d$, where $H(z)=\\dot{a}/a$ and $a=1/(1+z)$ is the expansion factor normalized to unity at $z=0$. \n",
    "\n",
    "Thus, the total distance traveled by photons to accumulate redshift $z=\\int_0^z dz$ is\n",
    "\n",
    "$$\n",
    "d_{\\rm c}=\\int\\limits_0^z \\delta d=\\frac{c}{H_0}\\int\\limits_0^z\\frac{dz^\\prime}{E(z^\\prime)}, \n",
    "$$\n",
    "\n",
    "where $E(z)=H(z)/H_0$ is the dimensionless Hubble function:\n",
    "\n",
    "$$\n",
    "E(z)\\equiv H(z)/H_0=\\sqrt{\\Omega_{\\rm m,0}\\,(1+z)^3+\\Omega_k\\,(1+z)^2+\\Omega_{\\Lambda,0}}\n",
    "$$\n",
    "\n",
    "The ratio $c/H_0$ is often denoted as $D_{\\rm H}$ -- the Hubble distance:\n",
    "\n",
    "$$\n",
    "D_{\\rm H}\\equiv\\frac{c}{H_0}\n",
    "= 3000\\,h^{-1}~{\\rm Mpc}= 9.26\\times10^{27}\\,h^{-1}~{\\rm cm}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Physical vs comoving distance\n",
    "\n",
    "In cosmology, the actual physical separation between two objects at a given redshift is called  *proper* or *physical* distance, while distance between these objects at $z=0$, assuming they were moving only with the Hubble expansion flow, is called *comoving* distance. \n",
    "\n",
    "Distance $d_{\\rm c}$, as defined above, is clearly the distance between observed object and us today ($z=0$) along the line-of-sight to the object because it computes the distance travelled by photon between an object and observer. Thus, $d_{\\rm c}$ is called the *line-of-sight comoving distance.* It is often also denoted by $\\chi$ in cosmology textbooks.  \n",
    "\n",
    "Photon's wavelength is increased (redshifted) by a factor $1+z$ between emission and arrival. This means that distances between all freely expanding objects in the universe increased by $1+z$, so relation between physical and comoving distances is given by\n",
    "\n",
    "$$d_{\\rm com} = d_{\\rm phys}\\, (1+z).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Angular diameter distance\n",
    "\n",
    "In astronomy we often want to evaluate physical size of an object or of a region that has redshift $z$ and subtends an angle on the sky $\\theta$. \n",
    "\n",
    "The *angular diameter distance* $d_{\\rm A}$ is defined as the\n",
    "ratio of an object's *physical* transverse size, $s$, to its angular size $\\theta$ (in\n",
    "radians): $s\\equiv \\theta\\, d_{\\rm A}$, where we assume small-angle regime. By definition, *we should use the angular diameter distance $d_{\\rm A}$ whenever we need to compute size of a region from its apparent angular size.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "How does $d_{\\rm A}$ relate to the line of sight comoving distance introduced above?\n",
    "It's easiest and conceptually clearest to work this out in the case of geometrically flat cosmological models ($\\Omega_{\\rm k}=0$), when usual Eucledian geometry applies.\n",
    "\n",
    "Let us consider 3 photons emitted at the same time from the center of the region, and from the distances $\\pm s/2$ from the center. As these photons travel towards us, the distance between them grows due to expansion of the universe, so that by the time they are observed by us the distance between off-center photons appears to be $s/a$, while the central photon travelled distance $d_{\\rm c}$ to us, as discussed above. Thus, the angle subtended by the emitting region will be \n",
    "\n",
    "$$\\theta=\\frac{s}{a d_{\\rm c}}\\equiv \\frac{s}{d_{\\rm A}}$$, where \n",
    "\n",
    "$$d_{\\rm A}=ad_{\\rm c}=d_{\\rm c}/(1+z)$$\n",
    "\n",
    "is an effective distance for which the usual small-angle formula applies. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In the case of non-flat cosmologies, relation is more complicated because space is curved, so we cannot use simple Eucledian expression relating angle and distance. In the general case, the angular diameter distance is given by: \n",
    "\n",
    "$$\n",
    "d_{\\rm A} = \\frac{d_{\\rm M}}{1+z}, \n",
    "$$\n",
    "where \n",
    "$$\n",
    "d_{\\rm M} = \\left\\{\n",
    "\\begin{array}{ll}\n",
    "D_{\\rm H}\\,\\frac{1}{\\sqrt{\\Omega_k}}\\,\\sinh\\left[\\sqrt{\\Omega_k}\\,d_{\\rm c}/D_{\\rm H}\\right] & {\\rm for}~\\Omega_k>0 \\\\\n",
    "d_{\\rm c} & {\\rm for}~\\Omega_k=0 \\\\\n",
    "D_{\\rm H}\\,\\frac{1}{\\sqrt{|\\Omega_k|}}\\,\\sin\\left[\\sqrt{|\\Omega_k|}\\,d_{\\rm c}/D_{\\rm H}\\right] & {\\rm for}~\\Omega_k<0\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Luminosity distance\n",
    "\n",
    "Another commonly encountered situation is when we need to estimate intrinsic luminosity of an object, $L$, given its observed flux, $f$. Luminosity is energy emitted per unit time, while flux is radiation energy received by a detector per unit area per unit time. \n",
    "\n",
    "In the non-expanding space, these are related by the distance as $L=4\\pi d^2 f$. In the expanding space, expansion affects the rate at which photons are received and the photon energy via redshift. \n",
    "The *luminosity distance* $d_{\\rm L}$ is defined so that $L=4\\pi d^2_{\\rm L} f$ in the expanding universe. By definition, *we should use the luminosity distance $d_{\\rm L}$ whenever we need to compute intrinsic luminosity of an object from its observed flux.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For simplicity, we will consider luminosity of photons of a given wavelength $\\lambda$; luminosity of photons of other waveleghts can be treated in the same way. \n",
    "\n",
    "Distance between two photons emitted at times separated by $\\Delta t_{\\rm e}$ will initially be $c\\Delta t_{\\rm e}$, but will grow as the universe expands by $a$ by the time the photons are observed. Thus, the photons will be received with larger time interval, $\\Delta t_{\\rm o}=\\Delta t_{\\rm e}/a$, i.e. their flux is reduced by a factor $a$. In addition, the wavelength of the photons, and thus their energy, decreases as they propagate to us by $a$. Thus, the flux through the sphere of radius $d_{\\rm c}$ -- the distance photons travelled to us - is reduced by $a^2$ compared to no expansion case: \n",
    "\n",
    "$$f=\\frac{La^2}{4\\pi\\,d_{\\rm c}^2}\\equiv \\frac{L}{4\\pi\\,d_{\\rm L}^2},$$ \n",
    "\n",
    "where \n",
    "\n",
    "$$d_{\\rm L}=\\frac{d_{\\rm c}}{a}=d_{\\rm c}(1+z)$$ is the *luminosity distance*.\n",
    "\n",
    "For a general (not necessarily flat) cosmogical model it is given by \n",
    "\n",
    "$$\n",
    "d_{\\rm L} = (1+z)\\,d_{\\rm M} = (1+z)^2\\,d_{\\rm A}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from codes.cosmology import d_l, d_a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "z = 0.10; Om0 = 0.30; Oml = 0.70:\n",
      "d_L = 3.222e+02 /h Mpc; d_A = 2.663e+02 /h Mpc\n"
     ]
    }
   ],
   "source": [
    "z = 0.1; Om0 = 0.3; OmL = 0.7\n",
    "d_H = 2997.92\n",
    "print(\"z = %.2f; Om0 = %.2f; Oml = %.2f:\"%(z, Om0, OmL))\n",
    "print(\"d_L = %.3e /h Mpc; d_A = %.3e /h Mpc\"%(d_H*d_l(z, Om0, OmL), d_H*d_a(z, Om0, OmL)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "zmin = 0.01; zmax = 10.; nz = 100\n",
    "z = np.linspace(zmin, zmax, nz)\n",
    "\n",
    "dL = np.zeros_like(z); dA = np.zeros_like(z)\n",
    "for i, zd in enumerate(z):\n",
    "    dL[i] = d_l(zd, Om0, OmL)\n",
    "    dA[i] = dL[i]/(1.+zd)**2\n",
    "    \n",
    "# convert to /h Mpc\n",
    "dL *= d_H; dA *= d_H"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is how to use the distance routine in the vectorized way. In this case, for many points the routine computing comoving los distance constructs cubic spline and then uses this spline to compute distances quickly for all points. This will be handy and efficient when we will need to evaluate distances for hundreds of thousands of galaxies in the SDSS. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# when called with a vector z of size larger than ninter, routine computes distances \n",
    "# using cubic spline constructed over [zmin, zmax] interval of z\n",
    "\n",
    "dLv = d_l(z, Om0, OmL, ninter=20)\n",
    "dAv = dLv/(1.+z)**2\n",
    "    \n",
    "# convert to /h Mpc\n",
    "dLv *= d_H; dAv *= d_H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlIAAAI0CAYAAADWcsc0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAimgAAIpoBvt37KgAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl8XFd99/HP0Uij2SSNNkvWYsvyHieOl4SExJQkTcJaCsEJT9m3xEChwANZKG1ZW0gotAVKSVKe0pBCIQ47AWIHErInjuMsTrzL+6J1NDOafeY8f9yRZHmRZWsZyfq+X695jebOuXeOHVn65pxzf8dYaxERERGR01dU6A6IiIiITFUKUiIiIiJnSEFKRERE5AwpSImIiIicIQUpERERkTOkICUiIiJyhhSkRERERM6QgpSIiIjIGVKQEhERETlDClIiIiIiZ0hBSkREROQMKUiJiIiInCEFKREREZEzpCAlIiIicoYUpERERETOkIKUiIiIyBlSkBIRERE5QwpSk4wxJljoPoiIiMjIFBe6A9OVMeZWYDUQAnYBQaA1/3plAbsmIiIiI6QRqcJrxQlUVcCt1lqFKBERkSlCI1KFdZW1dlehOyEiIiJnRiNSIiIiImdII1JjxBizArgVZ4quO/98u7X2jmFOW22MuTDfnvw511trQ+PaWRERERkTClLDyC8If9pau/YU7W4AbseZqlufPxYEHjDGXGutveokp8611l571HVuAtqMMXMUpkRERCY/Y60tdB8mFWNMK3AlcDPOQvCbrbW3DdN+BfAMsObY0ad8mOoBbrPW3jzCz7fA2qMDloiIiExOWiOVZ4xZbYzpwRlZCuKUJBiJOwFONIWXH1VaD9yUD2gjEcIJciIiIjLJKUjlWWvXWmsrrbVX5UegTjm1lh9xWgFsHKbZuvzz6qPOu8kY06PimyIiIlObgtTo9I8cdQ/Tpn9k621HHas+0Tn5YBUEfjImvRMREZFxpcXmo9O/iHy4acD+ka0VRx37Cs76q2PdA+yy1q4Zg76JiIjIOFOQGp3+MDTcNODAyJMxJmitDVlrQ8aY64FbjTFV+beDwLrhFraLiIjI5KIgNTpVp25yXPsQDCxE18iTiIjIFKYgNTojWSx+9GhVQRaXG2OeBRYCLiA+glMO5x8iIiKTyTwgALRZa5cXujOgIDVdzAG8+a/dI2hfgRO8REREJqM5he5APwWp8Xf0KFShqpVHgQqXy0UgEDhl4/r6eurr68e/VyIiIqdh06ZN9Pb2gvN7bVJQkBqd0w1Gw5VJGE87gMZVq1bx4IMPFqgLIiIio3PZZZfx0EMPgfN7bVJQHanR6Q9GI1r7pP3zREREzi4akRqdjThFOYfb/qX/vZFuOTNuHnnkEYLBoZkvFFK2ExEROVMakRqdp/PPw5VB6H9v/Tj3RURERCaYgtQoWGvX5r9cMUyzlfnndcO0mRCrVq0iFAoNeYiIiMiZU5A6tepTvH8HgDHmypO8fx3Oti9rT/K+iIiITFEKUifXeszzCeX3xdsF3H7se8aYm3AWoquCuYiIyFlIi83zjDErcDYNrmLoXXirjTEWp9RBN7DyBHffrQTuMcbsBG7Ot70WZzTqKmvtpFgfpcXmIiIiY0tBKs9au9EYc6KQBDgbDufbHfd+/thV+TB2Ac4o1j350SoRERE5SylIHWW4Ok8jqQFlrd2IUxJhUlJBThERkbGlICUiIiKTgrWW1JEUyT1JEnsSA4/kniTJw0nCL0cK3cXjKEiJiIjIhLA5S+pQisTuxNDHngTJjhRZoHhGKSXNpZQ0llLSVErxuQG8TaUU17jxvDUAjxf6TzGUgpSIiIiMCZuzpA6nSLQljgtLqa60E5Tq3ZQ0e5yQ1FyK+6IK/E2luCqmZiSZmr2WM6K79kREZDSstaS70k5Q6g9LbQnibXGSHWlyOUtxfakTlGY5z6UXV+BvLsVVNgaRw47+EmNNQUpEREQGZOPZwYC0K05iV/7r/QmyaYur2u2EpFkeZwpuWRm+Zg+u4OgjRS6eJX0gSWZ/kvS+JOn9ziNzIAFJS7ItOgZ/wrGlIDWN6K49ERHpX9Cd2DkYlOK74iTaEqTCGSh1UTK7lJLZHkpmeXBfVknFuzxU17tH/dm5WNYJSPsSTkDan//6YBKTtbjLi/G0ePDM9uCf7aF0hR/PbA+lzaW4PC7KLiuDh8bgL2EMKUiJiIicZXLpHIndCeI7405g2hknvjNO8kCSTMbiqivFPdvjBKZWL97LKylv9lDkGd2GJ7l41hlF2pvIByYnKGUOpiiyFncwH5RaPJTN9uC5pBxPiwd3o5ui4qm52YqClIiIyBSUjWWdkaSdCeI74s5jZ5xkR4qcyzijSbM9zsjSnwUJvmsmxaMcVbIZS/rAUUFpb8IZUdqXpCiTo6S8GM8cD945XspaPHguzQelhqkblE5FQUpERGSSysaygyFpe3zg62RPGjwuSubkw9IcL+431+Jr8eAqH92v9kxXmvTufEDakw9MexLYvgzFniK8LV48czxUzPHg+YsqvHO8lDaXUuQ+O4PSqShIiYiIFFA2kSWxM0Fse4z4NicwxbbHSIUy4HXhbvU6ganFg3tlOYE5Hop8rjP+PJvKkd6XJLUn4Ywo7XaeM4dTuFxQWu/G2+rF1+rBe0k5njkePHM8FAcUGU5EfyvTiMofiIgURi6TI7knSWybE5Zi22LEt8dJdKTAXeSMKM3x4J7jxf220YelbDjjjCbtTjiBabczumQjGUoCLmf6rdVL2VwvnlXleOd6cde7MUVmDP/U04OClIiIyBiw1pLuSBPbGiO2NR+YtsZI7E+StVDS4sHd6nVGmN5YQ3Wrd1TTcJmuNOm2OKnd+aC0O0Fqb5yiLJRWl+CZ68E710vlygDe62rxtHooHotaTjKE/kanEZU/EBEZvVwyR3yHE5JiW5zQFNseI92XwzXDjXuu13lcXEHw7fUU1535Au9MZ4p0W4JUW4JUW5x0m7N2qQiLp74U73wvgXlevFcG8c714mnxTNu1SoWiICUiInIC6a40sS0x+l7ucwLTFmd0KVdSNBiW5nnxrgpSMduDOcMAk+1ODwlKqbY46b1OWPI2ePDO81I234v39VV45+cXdp+ld8BNRQpSIiIybVlrSe5LEns5H5hedgJTqjdDUVUJ7nk+3PO8uC8JUvnumRTXntnoUi6eJbUrTmpnwpmO2xkntTtOUeaokaUFXnyvqcQ7v4HSWQpLU4WClIiInPVs1hJvixN7KUbfS33EXooR2x4nHc9R3FhK6QIv7vk+3KtnEJjrpch7+gu9bc46d8PtjJPelQ9LO+PYSBZ3pQvffB/eBV6qX1mO9711eOd4KSpVWJrqFKREROSsYbOW+M44fZv7BkJT37YYmYzF3eKEpdIFXrwfDFLR4sGUnH6QyUYypHbEncdO5zlzIEGxuwhvqwffQh/lC3348lNxxaOs6ySTm/7riojIlGNzlsSeBH0v9jmPzX30bY2RzVhK5jhhyb3Qh/+ySoKzPRjX6d/Wnz6YJLU9NhCaktvj0JehtKYE3yIfgUU+fG+swrfQR2lTqUoHTFMKUiIiMmlZa0kdTjlh6YV8aHq5j3Q8R8lsJyyVLvThf1UlwZbTD0w2Y0ntjpPaFie1I0Zqe5xUWxyXAe9sD77FPqoW+fC/rgrvAq/KB8hx9B0xjaggp4hMZploZjAwPd9H9MUoye4MrvpSShf6KF3sw/POesrnek/7DrlcMueMLG13wlJye4zMvgQlXhe++V585/iourgc3/vq8czxUHQGU34yPSlIiYjIhLK5/Dqm5/uIPh+l7/k+Ym1xrNdF6Tl+Shf5KL26ihl/00xR4PQWfR8dmJJbY6S2xsgcTOIOuPAt9uE/x0/N1UH8n2jUdJyMCQWpaUQFOUVkomUiGScwPRd1Hi9GSfflKGn1UrrYT+liH+VvqqW6ofS0rmvTOVK7EgNhKbk1RuZAArc/H5iW+Kl9fRW+TzdR2liKMQpMMj4UpEREZNSstST3J4luig48YjvjWJ8LzxI/pUv8lL51BvWfaTntW/5TexOktsRIbukjuSVGui1BscfgX+TDf66fmquC+D+ZH2FSYJIJpiAlIiKnxWYtsW0xos9GiT4bJfJshERHmuKGUkrPdUJT+WurqW7ynNZ1sz1pkltiJF/uI/myM8rkylp887z4z/MTvKgM/wfq8c7xntFdeCLjQUFKREROKpfK0be5j+jGKJGNEaLPRUn15XDP81J6bgDPxRXUfKABV8XIf53YjCW1I0bypfwo0+Y+cp1pPDPd+M/zU7U0gP/1VfgW+3B5Tr8wpshEUpASERHAWagdfT7qhKZnIkSej5JJW0oX+/GcF6D09TXUfWrWaVX9zvakSWzuc0LTS32ktscoLjYElvjxL/Uz4801+P9hNu4z3HpFpNAUpEREpqFcKkffi31ENkScx6YomRzOeqalATyrZ9Dw2ZbTKjOQ2h0nudkZYUps7iN7JIW3zk1gWYCa8wP431qNb6FPpQXkrKIgJSJylrNZS9/LfUSeHgxN6UTOGWk6Px+a/q5lxNul2IwltS1G4sU+ki9GSWzuwySy+Of7CCwPUP2aSgI3N1Naf3p34olMRQpSIiJnEWstyb1Jwk+FCT8ZJrIhQjKUwb3Qj2dZAM8baqi/cfaI75zLJXLO4u8X+ki8ECW5JUaxCwLn+gksDzDzvfX4l/pV8VumLX3nTyOqbC5y9sn0Zgg/HSb8RJjIUxH6dscpbvbiWR7Ae3GQGWuaRlzUMpfIOdNyL0RJPB8luaUPt7eIsmVlVCwP0PR/m/Av8Z92+QKRs5mClIjIFGGzlr7NfYSfyAen56LkvC68ywN4lpcR/EItNfUjW7RtUzkSL/WR2BQl+UKU5MsxSnyGsuVlVK4so+wNs/At1nomkVNRkJpGVNlcZGpJd6cJPxGm97Fewk+ESXSmcS/y411Rhuft9ZR/0Teieko266xpim+KOsHpxSglpUWULQtQeUEZZTfPwneOQpPImVCQEhGZBGzOEtsaI/xYmN5HewlvikCgGO8F5XhWllHzjpm4giP7kZ3elyD+bJTEsxESz0UpSucoOy9A+YVlNH6sgcDSgKbnRMaIgpSISAFkE1kiGyL0PtJL7yO9xPYmKZnnw3thGd6319H0xdYRjTZlwxkSm6IkNkaIPxshdySFf76X8leUM+P/1FL2tVYtBBcZR/rXJSIyAdI9aXof7aX34V56H+slFcvhWV6G98JyKr80l9oRrG2yOUtqa4z4MxHiGyOkXu7DU1VC2SvKqLu4nPKPN1J6mpv/isjoKEiJiIyD5KEkvQ/3EvpTiPCTYXIeF94Ly/FeVE7d9Y0U+U99J102lCG+MULimQjxDWGIZig/v4yKi8tpvmUWviU+ioo1RSdSSApSIiJjILE3QeihEKGHQoQ3RjFVJXgvLsf3mmoabpo9omKXqZ1x4hvCxJ+OkNwcxVNTQvnF5TRcGaT872ZRUlUyAX8SETkdClIiImcgsTdB6I8hQg+GCG+MUNRQivfiCnxvq6PpC62YouHXN9lUjsTzUeJPR4g/FSZ7MEnZuX7KLymn8dNN+M/za7RJZApQkBIRGYHkwSShP4bo+WMP4Q0RimaW4n1lBb531NH05VMHp1w0S3xDmNiTYeJPhSlK56i4qJyaVRVUXF+Pp9kzQX8SERlLClIiIieQ7k47wekPPfQ+HobKEnyrKvBdN7IRp0x3mviTYeKP9xLfGKE04KLi0goaX1NJxRdaKK7Qj1+Rs4H+JYuIANl4lt5He+lZ30PPgyGyLoN/VRDfa6tpvKUFU3yK4NSRIva4E5wSmyJ469xUvKqCGe+qo+xb83B5RrZNi4hMLQpSIjIt2Zwl+lyUnnU99KzvId6ZxntxBf5XVTDzhkaKvMMHn0xnithjYeKPOcHJ11BK8LIg9R+aSdnKBaoSLjJNKEhNI9q0WKa75OEkPff30P37bsLPRXGfG8D/Z0Gqvj6P4prh6zhlQxlij/USe6yXxDNhvPVugpcFmfmRBgIrA1oYLjJNKUiJyFkrl84RfixM12+76PljiJzXhf+yIP7rG6lc6Bv+3HiW+BNhYo/2EnuiF0+wmODlQereX0/Zd+drxElEAAWpaUWbFst0kNifoPu+brp/1010ZxzvK8rxX155yuk6m7Ukno3Q93AvsUdCFAOVrw7S9JZqKv65FdcppvpEZHpSkBKRKS2XyRF+Ikz3b7rpfqAbW16C/88rKbtpNtWzhy8pkNodJ/ZQiL4/hcgeTFL5ygrqr6okeEuzil+KyIgoSInIlJPuSdP9+266ft1F+LkongvL8V9RyczrG4YddcpFs/Q9GiL2YIj4M2HKFviourqKlm/Px9vqncA/gYicLRSkRGRKiO+M0/nLTrp+1UUylsV/RRWBDzRQucg/7HmJzX30/bGHvgd7KMZSdVUVc26YSdkrFmiBuIiMmoKUiExKNmeJbIjQ+fNOuu7vxtS6CbymiupvLqC4+uTTbtlIxpmu+2MPyReiVKwso/a1VSz8v02UDHOeiMiZUJASkUkjl8oReihE58866X44ROnSMgKvraLh+kaKPCcfPUpuj9H3QA/RB3ooMVB1dSUzP9VMYHnglBXIRURGQ0FKRAoqG8/Sc38PHfd20PtsFO+qCgKvrWbWZ1tOGoJsOkfsiTB967qJP9lL+ZIA1W+oYv7HluA+RT0oEZGxpCAlIhMu25el6zdddKztILIthv/ySgLvrid4a+Dk54Qy9P2hh+j93WT3xqm6opJZ75hBxXdU00lECkdBSkQmRCaaofs33bTf0050ZwL/VVWUfayJqnknL4yZPpAk+rsuor/vprjIUvOGGppunYP/FAvMRUQmioKUiIybbDxL933dtP+4ncj2OIGrqyj7xCyq55681EBya4zIfV30PdCNb6ab2rfUMu/eJbhnaMpORCYfBSkRGVO5dI6edT0c+eERwi/24b+yirKPNlE1/+QjT4nnokTu6yL2YA/lS/zUXVND1d/PorhMP6JEZHLTTykRGTWbs/Q+2kv7D9vpfjiE79WVlL13JpVLTr7mKb4xQuTXncQf6SV4QRkNb62h8rY5uDzaikVEpg4FKRE5Y30v93Hk7iN0/KoT97Jyyt9Sw+x/mHPS9vGNESK/7CT+aIjgxeU0XzuD4DfnabG4iExZClIiclpSHSnaf9jO4f89AtVuyt9aS/Ovzse4TxyGEi9Eifyyk9hDIYIXldF83QyC31Z4EpGzg4KUiJxSLpWj6zddHP7vw8QOpym/ppb6/zoHV8WJf4SkdsYJ/7yDvnXdVJwfoPFtM6j8xlyKThK2RESmKgUpETmp6AtRDv+/w3Su78F3RSUVn55NzUnuuMscSRH+eQeRX3USmO2h7u0zqP5SC65hNhEWEZnqFKREZIhMb4YjPzrC4buOQE0JFX9Vx6xPzzphlfFcLEvkvi7C93bgdkHd2+tY8MD5lAS1p52ITA8KUtPII488QjAYHHIsFAoVqDcymVhrCT8R5tCdhwhtjFL+1lrq7lyEq+rEgajv4RDhe9rJ7ooz49paWu5ahGeWZ4J7LSJSeApSItNYpjfD4bsOc/iuIxTN9RJ8Rx2z/3HuCdum2uL0/qSd2Ppuqq+oZP4tswisCGCMNgUWkelLQWoaWbVqFQ8++GChuyGTQGRjhIP/cZCep8OUX1dH/d3n4Co//sdBLp4l8qsuen9yBF9VCTPfV0/1ra0UlWrRuIgIKEiJTBu5ZI72H7dz4PaDMMNN8N31zP5i6wnbxjdF6P3hEdKbo9S/rY459yyhtKF0gnssIjL5KUiJnOWSB5Ic+I8DtP+8E//ra5jxHwsprj1+37psJEP43g7C97RTca6fOTfMpGJVhabuRESGoSAlcpbqfaKX/f+yn+i+JMH3zGTWb5dhio8PRYnno4TuOkxmax8z313P/PXnU3KSReYiIjKUgpTIWSSXydF5byf7vrUfU19K5Q0NVJ13/H53uXiWyC876f3hEcrme5nz4QYqLtXok4jI6VKQEjkLZCIZDt15iIPfP4z38kpmfGshxfXHT9+l9yUI/fdh4g/3MPMd9bT++jzcM45vJyIiI6MgJTKFJQ8lOfDNAxz5dScVb6+n6efnUXSCSuJ9D4cI/dchSrKWxr9upPqb8zAujT6JiIyWgpTIFBTbHmPfbfvo2RSl8oMNzP7tsuOCUS6RI/yzDnp/cJiqV5az+Otz8Z/jL1CPRUTOTgpSIlNI9Pkoe7+yl+iBJFUfbTph+YJMR4rQfx+m7/4uGt5Zz7wHtHhcRGS8KEiJTAHhp8Ls+fIekmlL1UebqFxRdlyb1M443d89QHZ7jOaPN1H71TkUlahwpojIeFKQEpnEep/oZc8X95AuKaL6k83ULDp+ai6+IUz3dw5QgqXl07MIXhHU3XciIhNEQUpkEup9vJfdX9hNxuei+ubZlC7wHdem78Eeuv59P2UtHhZ9tZXA+ceXORARkfGlICUyiUSejdD2D22kXUVUf7aF0nnHB6jIbzrpvv0AVReWs/QHi/HO8RagpyIiAgpSIpNC3+Y+2j7XRiKao+bGWZQec3edzVkiv+qk5/YD1F5ZxfJfnkfpTO19JyJSaApSIgWU2Jdg9+d2E94Vp+am2VQvH7qI3GYtkV900nPnQWa8rooV687HXaMCmiIik4WClEgBpHvS7P3KXjr+0EPNjbNp/rPgcW3Cv+ig5z8OMOON1ax84HxKqlXCQERkslGQEplAuWSOA98+wIHvH6byw43M/tWs49pEft9F9zf3U3tlJSse0AiUiMhkpiAlMgGstXT+rJO2L+7G/5e1zPr1Uox7aI2nvkdCdH1tL1UXlbP8N1oDJSIyFShITSLGmGeAm6216wvdFxk7kY0Rdnx6B6bFy8z/WULxMVXGE5ujdPzjHspmlXL+PUvwzPYUpqMiInLaFKQmCWPMPcAK4PjFMjIlpTpTtP1tG7074sz4fCul84eWMkjvT9D51b0UJ7Oc82/zCCxVHSgRkalGQWoSMMbcCuwqdD9kbOQyOQ7dfoh9tx+k+lOzaP7C0P3wctEsXd/cR2pDmLn/1Erl5ZUF6qmIiIyWNuIqMGPMaqAL+HGh+yKj1/tYLxsv3Uj3niTNv1hK4OqqgfdszhL6n8Psf+vzzLigjJWPrlCIEhGZ4jQiNUaMMSuAW4EqoDv/fLu19o5TnHOVtXZN/muZotI9aXbdsotwW4K6by+kpGnoOqfYE710fGk3dW+qZsHjK3H5XQXqqYiIjCUFqWHkp9yettauPUW7G4DbcULR+vyxIPCAMeZaa+1VJzgnCHzGWnvtOHRdJoi1lvYftdP2lT1UfWIWTZ8fOo2XPpSk40u78ZQYlv38XDyztJBcRORsoiB1DGNMK3AlcDPQmn8erv0KnBC15ui77ay1IWPMnwM9xphbrbXHXudO4Pox7bxMqPjuONs+tA2aPTT/dClFgcFRJpvK0X3nQeK/7WLe1+ZSeYWm8EREzkZaI5VnjFltjOnBCUVBRr74+06AE03hWWtDwHrgpnxA6/+s24Gv5N+XKcbmLAf+/QDPvflFyj7WzIwvtg4JUbEnetn7l89TUVXMyidWKESJiJzFNCKVl5++G5jCM8ZceKpz8tNzK4CNwzRbhzPCtRq4LX/OBUCVMebodv1B6zPGmLcBPz7VlKJMvNi2GFvXbKV4WRmzfjG0qGa2O037l3fjTuZY9qvz8DRrGk9E5GynIDU6V+afu4dp0z+y9Tbgtvwo1MpjG+Xv3rsHZ6RKAWqSsTnLgW8d4MBdh5lx61w8S4bWfOpd207vnQeZ++U51LyppkC9FBGRiaYgNTr9i8iHmwbsn7471V15Vcc8yySR2Jtgywe2ULQkQPPPzsOUDI5CpQ8mOXLzToJLfKx8bDnFAf2TEhGZTvRTf3T6p+OGW+s0MFpljAkeuy7KGNO/sP2C/KFbjTHXArdqq5jCstZy5O4j7P76PmZ8eS7elWVD3g/ddZjIjw+z8FsLqLi0okC9FBGRQlKQGp3THT2q4pjQlQ9LExKYNm3axGWXXXbKdu9973t573vfO+79mczSoTTb1mwj7XHRfO95FPkGF5On9yc4/Kkd1FxawYLHVuDyqCaUiMh0pSA1OiPZF+/o4FTQffR6e3t56KGHTtluJGHrbNb7RC9bP7SNyk/OourqoVm598dHCN91iEXfXUj5K8oL1EMREZksFKSmkYqKCpYtW3bKdi0tLePfmUnI5ix7b9vL4V93MfN7iylpKB14L9OV5siNO6iY52XlIytw+TQKJSIiClIT4ehRqILWjVq2bBkPPvhgIbswaaXaU7z87pcpOjdA8/+ei3ENlqaI/qGH7tt2M/+f51F1le4FEBGRQZMiSBljynHWD7XiBI8QziLtkLV2dwG7diqnG4yGK5MgBdL7RC8v37CV2n9oxXfUonGbytHxlT242lOs+ONySqpKCthLERGZjAoSpIwxy3DqKl3JKcoC5ItWbsQpbPkTa+2mce/gyPUHoxGtfSp0JfNHHnmEYHBoV0Oh6Vtc3VqnQvmB/2mn4fvnUDJzcCov1Rbn0Me30fy+mTR8eD7HFE8VEREBJjhIGWM+yOAedgan/tID+ecQ0JV/DgLV+edWnNIAK4GbjTG7gK9aa783kX0/iY04YbB1mDb97410yxmZANm+LFvXbCXlcdH84yVDakOFf9ZB+HsHOfe/FhFYGhjmKiIiMt1NSJAyxlwB3IETKtYDt+GMLvWexjUqcApgrgHuNMbcAtxgrf3jOHR5pJ7OPw+3cKb/vYLXhFq1apXWSOEU2Nx83WYC755J3VtqB47bdI72z7dRmras+NNyXH4tKBcRkeGN+6bFxpiv4oSItdbaImvt1dbaO08nRAFYa3uttWuttVcBlcBPgXuNMf80Dt0eaZ/6t3IZbnqyfzuYdePcHRmB3sd62fSG56n6QivlR4Wo9KEk+657kdoLylh01yKFKBERGZFxDVLGmPuBOUCltfaWsbpuPlT1TxHOM8b8fqyufQLVp3j/DhioUH4i1wG7tH9e4R3+78Ns+eQOGn+wBM9RU3axR3s59L6XWfSv82jfQT6xAAAgAElEQVT8SKPWQ4mIyIiN29SeMeYnwDpr7dfG6zPyi7evM8bcaIz5irX2M2N4+dZjnk/WhzX5EHU7MPfo94wxN+Gs87p2DPt1xqbrYnObs+z6zC5C2xM0/e+5FHkG//+h5z8Pknqoh+Xrz8dd6y5gL0VEZCoazzVSX7HWPjuO1x9grf2aMWb5aK5hjFkB3IOzpunotLHaGGMZLMmw8gR3360E7jHG7MRZTB/CCU/XAVdpz7zCySVzvPyel8k1e2j49wUDx206x5G/3YW/wsWi3y6lqGTcZ7lFROQsNG5BaqJC1Fh9nrV2ozHmRCEJcDYczrc77v38savyYewCnFGse6y1a0bTp7E23Rabp3vSvLj6Rbx/UUvV2+oGjme60xxcs4Wmt9fR+JHGAvZQRESmuklRkHOyGK7O00hqQFlrN+KURJACS+xL8MI1L1L5yVn4L68cOJ7cFuPwR7ex8N/mUXlF5TBXEBERObWCz2cYY64xxnzFGNNygvcqjDEb8vWnREYk+mKU5974AjX/OHdIiIo93kv7x7ex9KdLFKJERGRMFDRI5UsjrAVuwlmsPUT+7rwLgEpjzH9MdP9k6glvCPPiO16m/vZFeM4bvDMv/MtOQl/fy7L7z8c331fAHoqIyNmkYEEqX2DzSpzF2b04C71PKH/n3y6NTI1O/117Rz/OJqFHQrz0wa00/Ndi3LM9A8e7v3uA1G86WXb/+bhn6M48EREZO4UckboS+LG19mvW2ipr7X8O1zgfpibV4m2ZPLrXdbPl4ztovHvJkD3z2j/fRvHhJEt+ukRFNkVEZMwVMki1AjtP85yzawhlgq1atYpQKDTkcTbo/EUn2/++jaa7l1BcXQI4taMOf3oHZTXFLPjuAoqKC74cUEREzkKFvGsvxCmKXZ7AcHvayTTU+atOdt66l6a7l1CUH3GyGcuhj2+j5sIyZn9mdoF7KCIiZ7NC/m/6BuAzxpiykTQ2xnwlf44IAF2/62LnP+6h6b/PGQhRuWSOgzdsoe7yoEKUiIiMu4IFqXwBzd3AA8aY84drm7+77ybg1gnomkwBPQ/0sOPvd9P4/XMoCuRDVDzLgfe9TOO1tTR9rKnAPRQRkemg0AU5r8cZZdpojNmY/3onzrTfXJypvytx1kbdaa39Q6E6KpNH6OEQW2/aSdNd5+CqcL6Fc8kcBz64hVkfqKf+nfUF7qGIiEwXBQ1S+W1ZLsApfbASWHFME5N/vs1ae8uEdu4sdDZsWhx+OsyWj+2g8QeLcVXlF5ancxxcs5Xmd9YpRImIyIQq+K1M1tqN1tq5wIeAB3BqSpn881qcTYIVooTYjhgvvX8rDd9bRHGNUw/KZiwHP7KNhjdXM/P9MwvcQxERmW4KPbU3wFp7B3BHoftxNpvKmxanjqR4cfVm6r+1gJIGp06UzTl359VfVUnjh7X5sIiITLyCj0iJnEommuGFa16k+vNzKF04uL3LkVt2UvOKMpo+roXlIiJSGApSMqnl0jk2X7uZsvc34LuoYuB45zf2UlZTwqybZxWwdyIiMt1NmiBljPmgMebHxpjtxphs/vk7xpiWQvdNCsNay7YPbaPkzyope331wPHQ/xymaF+CuV+fizFmmCuIiIiMr4KvkTLGLMO5a6+Vwbv0YLD8wRpjzA3W2u8Von9SOPv/dT8JY6h77+Ai8ui6bpK/62Lpb5ZiihSiRESksAoapIwxc4BncAJUCFgP7AK6cILUBTglEe4wxgSttV8vVF9lYnXf382hX3TSdPeSgWPxZyL0fns/y9adT1HppBlMFRGRaazQI1L34JQ5uN5ae++JGuTD1m3AbcaYe621uyewf1IAse0xtt24k6YfLcG4nFGn9IEk7Z/ZwfLfn09xRaG/bUVERBwF+41kjFmOM3W3YrhwZK1tA641xtwD3Ax8eGJ6ePaZCgU5M70ZNv+fl6j/5gJclU7BzVw8y6E1W1jy/cWU5ksfiIiITAaFnB+5ALhjpCNM1tpr8+fIWcpmLS+98yWCf9NM6YLBMgeHP7mDlk83U7ZyRPtbi4iITJhCLzR56jTb94xLL6aJVatWEQqFhjwmkz1f2YNZ5CdwddXAsa5/20flOT7q/qqugD0TERE5sUIGqfXAhWNxofydfzKFhR4O0X5/D9WfaB44Fr2/G/tSH3O+PKeAPRMRETm5ggWp/NqnlcaY2SNpb4z5NM7i9BN5Zsw6JhMu1Zliy19vY+a3FgyUNEjtjBP65j4W3714YMG5iIjIZFPoqb1rgbXGmPNP1sAY02KM+QrO6NWPj3mvPL9oXaYom7Nsed8Wav62heIZzkbEuWSOQx/fxjl3Laa4XHfoiYjI5FXIu/ZuBL6af7lxhBWqV6uS9dll/7/shxYv/ldXDhzr+EIbsz/SiH+Jv4A9ExERObVC/u9+iKGVzEfDjtF1ZAKFnwxzcG0HzT85d+BY5L4u3Kkc9e+vL2DPRERERqaQQWoXzl14c6y14TO9iDGmFdg+Zr2SCZGNZdnyoa3MvHMxpjhfdHNfgtC397PiT8u0h56IiEwJhVwj1Q20jSZEAVhrdzF2I1syQdr+ro2yd86kpNEpsGnTOQ79zTYW/7+FFJdpXZSIiEwNBfuNZa191hjz52N0ubljdJ2z2mSpbN77aC/dz0ZpvnHwhs3Ob+yj8e11lC1X0U0REZk6xn1EyhhTfrL3rLW9Y3GtfCkFmQKy8SxbP7qN+n+eN3As8XyU3OYojR9tLGDPRERETt9ETO2NZchRYBqFyVDZvO3v2yh7Rz0l+T3zbCrHkc/sZOEdCwdqSImIiEwVExGkKk/dpCDXkgnW+1gv3RuiBN8xeEde57/so/kDM/HO8RawZyIiImdmQtZIGWM+BTw7ysusRGUOpqxsIj+ld8figWOJF6LkXozS8A0tcRMRkalpohab3zZBnyOT1P6v7ydwzYwhd+kduWUn569doik9ERGZsiaq/IEZo4dMQYm9CQ7d00HwPTMHjnX9236a3l+Pt1VTeiIiMnVNRJC6BVgL3GCtLTrTB3DBBPRVxsHOT++k9u9bBgpvpnbFST3ZS+NHdJeeiIhMbeM+tWetvQ3AGPPnxpifAF3APdbaP5zmdTYalbuecnr+0EMiaam6pGLgWPvn21j4r/M0pSciIlPehBXktNY+ADxgjKkArjPGfAjYCdxurd09wsusH6/+ydjLpXPsuGkndbcvGjgWXddNeYuH8gtOWl5MRERkypjwyub5Ipx3AncaY+YAH8o/r7PW/ucpzr16IvooY+Pgdw7ivbJqoGZULpmj6+t7WbH+/AL3TEREZGwUdFOzfEXyW+C4qb/brbWbCtk3GZ1Ue4r93ztI8y8HQ1PP7Qdoun4m7lp3AXsmIiIydibN7rAnmPr7W+Ap4I7RbmwsE6/tH9qo+r+zKCp17mdIH0ySeKCHhseWF7hnIiIiY2eiyh+MmLW211p7p7X2OuAB4FZjzO+NMW8pdN+muv5Ni49+jIf4rjihZ6OUvaZ64FjHl3Yz72tzKSqedN9yIiIiZ2zSjEgdyxjTAlyXfwSBK40xu6y18wvZLzm1PV/aQ/WnZg28jj8TwVNqCL56fIKbiIhIoUyqIGWMKccJTmuAFf2HgV3A7cAdBeraWWHVqlU8+OCD4/oZsa0xenfEmfWqwdDU+bU9nHvnwnH9XBERkUKYFEHKGHMFTnha3X8o/3wHzsLz0e7TJxNk9+d3U3vj4GhU7LFeyud68S3wFbBXIiIi46NgQSo/dbcGuAFn6q4/PG3ECU93FqZncqaiz0fpO5Ki8sLBGlFd39jLeXcvHuYsERGRqWvCg5Qx5hrgMwydugsxOPrUNtF9krGx+3O7qbl59sDrvj/2EFwawDtH++mJiMjZaUKClDFmGYOjTzA4+rQWJzw9MMLr/N5a+5px6KKMUnhDmEQiR/XSwMCxzn/dx7KfLilgr0RERMbXuAcpY8zTOKNPQ6bugJ/kq5yfjivHsm8ydnb/w25qbhpcGxX5XRc1l1bgafIUsFciIiLjayJGpFYCPQxO3e0+k4sYY747lp2SsRPZFCFVZChd5B841v2t/az47dIC9kpERGT8TdQaKQNcC6w2xpyq7YnObc1/bceyUzI29v/LfqrWNAy8Dv+yk9qrK3HXaSsYERE5u01UkLpitHvnGWNuBz44Rv2RMZI8nCT8coxZFwzeqRf63kGW//q8AvZKRERkYkzIfh1jtAHxzWNwDRljB79zkOD7Zw68jj0VpuJcP+4ZGo0SEZGz30QEqfVjcRFrbQhn7z2ZJLLxLEd+1knZG2oGjvXccZCmTzYVsFciIiITZ9yDlLX26sl4LRm9I/9zhLI312CKnXVvqb0JilNZAucFTnGmiIjI2WFCpvbk7GOt5cB3D1LxjvqBY6HvHaL5E80F7JWIiMjEUpCSM9KzroeSpQFcFc79CtlIhuSGXqpeV1XgnomIiEwcBSk5I/v/dT+VHxgsedD7o3Yar2/AFJ12eQsREZEpS0FKTlvfy30ks+BudfbQs1lLZG079e+uP8WZIiIiZxcFKTlth+44ROX7BkNT9LddzHhjNS6/q4C9EhERmXgKUtPII488QjAYHPI4Xbl0js713fheNXhuz/cP0fDXDcOcJSIicnZSkJLT0nN/D75XVw6shUrtjuOrKtHmxCIiMi0pSE0jq1atIhQKDXmcrsN3HaZ89YyB1+G1HdS9q24suykiIjJlKEjJiKVDaSI745Qu9A0ci/6+i+q/qC5gr0RERApHQUpGrGNtB+V/WTvwOvZkL1WXVODyapG5iIhMTwpSMmJH/ucIZW8ZDFLhezqoe7em9UREZPpSkJIRie+Kkykpori6BIBcIkfqhSgVl1YUuGciIiKFU9AgZYy5opCfLyN35O4jVBy1yDx6fzczrqnBGFUyFxGR6avQI1JrCvz5MgLWWtp/1oH/6sF99ML3tFP/LlUyFxGR6a3QQaq1wJ8vIxB+PEzJeWUUlTrfLpmOFK5UFu9cb4F7JiIiUljFY31BY8yPgZGWzF4x1p8vY+/wXYepuPaoReY/7aD+nRqNEhERGfMghTPKdRVw+tUeZdKxWUvP42Fmf25w8DDyiw4W/GFZAXslIiIyOYxHkLoJ2GmtveVUDY0xG8bh82UMhZ8M41lZNvA6uT1GYK6XksqSAvZKRERkchjzNVLW2jZg+Qib7xrrz5ex1XVfF/4/H1xk3re+h5o31RSwRyIiIpPHeC02f3aE7daN0+fLGOlZH8J3yWCtqOgD3VS9rmqYM0RERKaPcQlSI5nWy7e7czw+X8ZG8kCSXLmLIo/zbZLtzVCMxV3rLnDPREREJodClz+QSazrvi4CR0/rPRSi+jXaoFhERKSfgpScVPd93fivrBx43fdAN9VvUJASERHppyAlJ5RL5ojuTVDSUAqAzVlSL/YRWB4ocM9EREQmjwkNUsaYayby8+TMhf40dJF5YlOU4MXlmCLtrSciItJvPOpIDedO4KcT/JmTkjEmCHyGwSrwrThFTL9ird1YsI7lOdN6R62PeqCHek3riYiIDDHRQUrDGYOeAdZYa9f3HzDG3A48Y4xZWegw1fNwiMZPzRp43fdgD5Wfn13AHp1cJpOlsz1COJygL5oceMRiSVKpLKlUhlQyQzqVJZvNkbMWLPlnS5GrCFeR4erXn8e8hXUn/IxUMsOGp9rweErwet34/G58Pjc+fylen5sijdSJiExLEx2k7AR/3qRkjGnFGYFaA6w/6q11wA3A24CCBanY9hhFszwYlxMOMu0pPJXFFJdP9LeLI5XKAOB2Fx9zHNxu2NPWxUc/cBcH2xdjTI6ZtVsH2ry883KKigwL5zw+cOz5ra8DYOnC3w4ce3H71fzydx5+9/vB63/ucxCPw223QU9PH1/67C842L6IXK6YpvoXB9pFY5X4fMXU1yUIlHmcR8BDWfngo7zCS3mFlyXnNVIR9I3p34+IiBROYX4zTnPW2l3GmJUcX9n9wvzz0xPcpSG67+seWvbgDz1Uv35ipvWy2RxtOzrYtvUw27ceYcfWI2zaVMQ73/9q1nykaaDdvHkwfz789rdQVe0HoG3/BRS7UkOCVCjcgDG5IZ/RF6/k2MHRaF81e/eWDjm2fj0kk87XiXgagH2HlpHJuocEqZd2XEUq7WPVyu/TfiQCwHNb3kAm62blkp8NtDtw5ByueVsVn//iYJDq7ASfz3kcPBDil/dupLLKT7DSR2WVn8oqH1XVAYJBH65i3RsiIjLZKEiNEWPMCuBWoArozj/fbq2940Ttj526y59/A3CHtXbtOHd3WF33dVH9jfkDr6Pre2j+59Zhzhid9iNhNj69m41P7+aRP/XQ1VVKVXD/wPvPbP4gnd90seYjg+fU1cGhQ87XFUEfRS5DcXGKdNoz5NpFRRmy2aEFRA0Wa4cGKWuLcLmG9iuZBE/+cv1BKpstodiVHNIuk3VTXJwYciwWD4IZOgC779D5/PRn5Xz+i4PHVq+Gnh547jnYt6eLX6zdyPY9l2AIM2/24PBYX7yKYNBFY6OlqjpAVbU//8h/XROgOv/aXap/1iIiE0U/cYdhjLkVePpUwcYYcwNwO3BV/5qn/GLyB4wx11prrxrm3HtwFpxfCdxmrb15zP4AZyATzZCIZHFVOZsS21SO7J44vkVjOx3V3dXHH9e9xPrfvcTuXR0Dxx/e8D6MybFq5X8PHCsPHOHAgZYh5/ePRgEUFRkqK/1UB/eQSg3t5+zGjVg7dCRnQcvD2GNGpOa3PMx7rl8FNAwc+9jHoDj/LySRcIJUVXAvxcVDg5QxOUpLYkOOpTMe/L7uIcdSaR8z6oZ+bns7NDY6X4e6nWsc6ZxPSXGCebMHpyOf2/J6ikyWi5f9iN27OgFn2tLlSrKg5Q8D7cLRWv7fD/+KhYtOXH0+FIqRSWeprPLjcmmES0RktBSkjpFfv3QlcDPOOqZhg01+JOl2jlk4bq0NGWP+HOgxxtx6soBkrb02f50gcI8xZidOICvIhs6hB0P4VgUHXseeDFN1eSXGjH4xdSaT5fGHd7Duty/yzNO72bLzEnojr+bC8wZzak3lbvYcXE4mW0KxywkvleUHKHLVkEqV4s7ng69+1Xn0e9M1y3ntGzMEyiz+wGvxB0rx+Uop9RTjdhfjdl9CSYmL4hIXxhiMMRTlc0Qua8lmc5QHvUP6+573DH7dNKuKT9z8GhKJNPFYMfHYK4j1pYjFUlx0yZP0RZNEIzVEo0mi0QRLF/4GUzR0SrG2ahevfe3Q/bw7OmB5/lBPTx/gBK7AMSEsmfITLDs05Fh791wC3qHtnt/2F/z1R92sP2rl3Zo1MGsWfPaz8Mu1G/nRXU/QF6+httbQ3OyiuqaM6toANTUBamaUUV0TGHjt8Wo7IBGR4ShI5RljVuOUZ9iAs+h7F06QOpU7AU40hZcPU+uBm4wxtw8XjvJtrwV68p8/9/T/FKMXfjyM94KygdexP4VoHOUmxel0lvW/28wPv/80nR09Q97rCrUQjVUS8DnHayp309nTQirlo9jbi9tdzGuuOsSKC9243asGzquvH/oZ173zolH18VSqawK85g3njbh9Op0lEk4QCccJ98YJhxPE+pJc9bqho0AHDgyuw+ofkWptfgpvae9AG5ufISx19w05lkz5qa7YO+R6yVSAhoYhh/jZz+Dyy52vuzqjADz9/DWUlx1hRfQXwGEANu+4kjLfAWY1PDdwrimupWFmETMbvNTUllFTW0ZLaw2XvnrBiP8uRETOZgpSefnpu4GhEWPMhcM0728TBFYw/B1263BGuFYDt+XPW40z6rTmmD6EjDEhoNUY01qIUanwU2Fq3jNz4HX8mQjlX2wZ1TU/d/NP+e1vS9m84+28cvkP8HnCADTOeIld+17JwfZzWNDyKACtLYdZfe3jnL/iIhYurmdWSzXFxa7hLj8plZS4BtYxDcftZmCU7T3Xr+KN1ywj1B2jp6ePnq5ZhHpi9HT3cdGl6+ju6qOnO0BPdx/ZDCxq/SN+72AwzeVcuN1ZZg1WrSCVcka9+qcPuzqjWGtIpgJ43DsH2lkLB44soaH25SH9W/fw25hRtXPIHY5dfW/ixlvg2msH20UiEAiAMXD/fS+QSmXzwcsZ5aqo8I7JqKaIyGSjIDU6V+afu4dp0x+G3kY+SOGUPbgyP+V3bFgKgnNn35j1coRszpLsyeAqc74tbNZiYllK8uulztTlVy3m4T89RTrjYde+izh3/joA/L4eVi65lxUr4lx+5SpWvqKFeQvqpm1NJq/PTZOviqbm4UcAczlLuDdGd1cfXZ1RurvOoburj+7OKB/6m8285boLjmoL3/kOnJcfTOvqiJLOeHAVpyh1RwfapTMecrkSSksHj2WyJWQyniHtcrkiNjw7n0cfHQxSuRxUVztTiN/6Fqz94dPs3dPNjr2XUFWxgergPkrcLmfKsKaM2hmDAas/bNXWlhGs8k/b//YiMnUpSI1O/yLy4UJPKP+84qhjtwOhY8NSviAnwLUUQGxLDPeCwcXayS19lC09s731fvUr5xfsX/4lXHHVOfzorifYsXcLHV1zybb+keoaN5dftZgrX7uE1nkzxuqPMC0UFRmClX6Clf5T/t15PPDhDw++vuI157BgcTevvPTXdLZH6eryEI0kMCbHgpaHqKw4MNDWWbifw3NUuEqmAlhraBqsREFHB6TTg9OtnZ1RMlk3O/deAs2PUx3cRzqV5cC+Pu66573MnfU4c5ufAiBnDfsOnU9VxX6C5d3O2qzaMt68egWvunzhqP+uRETGmyqbj07/GqrQMG0GRquMMUFrbchau9YYsyt/xx44o1BVOIGsYFXNw0+G8SwfXB+V2Bil/OLy075OOg033ujUSLr4YqirK+Kv3n0xbW0PMmu2n/d84Gr+7IqFU3LKbqq79u2vOO5YIpGmqyNKV2eEzo5aujrn0dkRoasjyvnLfkRHRx/hkCGXs5iiLG99czuXXDIY4PbnK1U0NUFfX5J4LEUi6dQdKy2NDH5OKkAuV4KrKDNwLJkK8PLOK1nQ8hBl/k7aD4fZvy/N3Wvn8NWvwvve57SLx+E3v4FXvAKqKhP83afXDqzZ6g9fNbUBqmvLqKlRCQgRmTgT/dPmJxP8eePtdFdhV5EPXfmwNKEjT5s2beKyyy476ftv3PZGVv/ggwOv489GaPxU00nb97PWEutL4Q84BS1LSuD734dLL4VPfAJ+9CNnVMrrdfPKV83TbfeTjMdTQmNzJY3NlSdtk83k6O7uo7MjQl2dj6qawfeWL4cjR8Drhc4OJzgVFWVpmLGZcv9gaYtE0gnpnqPDVTJwgmNltLe7BxbZA+ze7Uwlfuc78NqrI2x9+TCbno2xacuFzG1+ghnVzuBuJltCODqDmfV9NDaWOEHrqLsRB9Zt1Zbh87u1bktERm1Cg5S19kMT+XkTIHjqJkNGq0bSftz09vby0EMPnfT9d/veg3vhUVN7L0QJnGJqL5vN8d1/+wM//F8/t31jOZeucipYXnyx80vv1a922rmKi1h1me70mqpcxUXUznDWNx2rqAhm5AeoSt1BvnXnu+jsiNDZEaGjvZmuzko62yMcOhjD572b4qLBBfJFJkd1cPeQRfP94ero6cOjR736w1o8WU5vpIFcbvDHWDRWzVPP/xXnJe4jGX+JXTs6iMYq2dp2Ga1NTw1MXWayJdz8+Xdz6apKvEOrXgDQF02STKSpqPQp+IvIsDT+PY1UVFSwbNmyE75Xki3hnOzSgf31suEMpRXFFLlP/ksknc7yz/94H+t/t5enn30/V7/GsuHpFIvPcW5DW7PmpKfKWcpdWsy8hXUn3fzZWksknBgIWp3tUTo79tPZUUNnh4fO9ghFrg5+9o8pLrjAfdR5sHQpzJ4NB/Y4a7b6A1fpMaNZMHSEK5aopKN7LrNmPjtwrDc8k6uuruTnP3fW8QG88AJ87Wvw6U/DvraX+fbX15OjhECgisYmFzNm+KmuDQyMbDmL552H16d6WyLTlYLU+Dt6FGq4tVTjbtmyZTz44IMnfC/0cIi9PzlqGua5KOWvOPn6qHgsxZf+7hc8u2EPbjesXPJTnn7hOl7/um627aihpETrn+R4xpiBDZyHWyhvreXoWberr3YeALu3ZwiUeahMHGTpwt8MKV5qc0WUlkSHBKnkwPTh4KL5rHX+WfaXhQDYvBl+8APnfwA62522HZ2zuf/hN3PheT+mOrgdgN7oDA4cPpeWxkfxeZ16X26Pj8qqADNnegbCVWNzJVe/fuS1x0RkaprUQcoYcyNgcfafCxe6PydwusFouDIJBXXcQvNnItRedPw0Djgh6sa/uZcXno/iy+9FFyw/xIolP6VhZoLOjr9kZkNBZzFlihtu7dKb3rqCN711xcAi+c6O8+hoj9DVGaWzPUJX53o62kvp6vQR6okRLD/IotY/4vUMFjkt9ZZRUjI0SB3I37DY2AhP/skJYolUftTrqBIQ4Ugdew+toHnm8wPH2va08Mt1r+eVy+6mouwlALxli3ngofP4xCdgZr40WyLhTIW63fD7X79AZ0fE2Sdx4OGnvMKnMhAiU8ikDlLW2q8ZY94K/MEY80/W2p8Wuk/H6A9GI0oN1tqCjkg98sgjBINDuxoKOV2KPBmh/LMtA8fjz0Yo/9BMjmWt5V+++nvu/cUKQuGZXHT+j/B6nF865y6J8E9fX029QpRMgJEskk+ns3R3RfOBa+ZA2PKXwf95F0NGvV75SmcbnYYG6OxwglMq7awZ9LiHloCAoZXmkycIXNFYE7fdNnSK+3vfg49+FLZuhT+se4nnn91Hd28jXT0tzG68D3dJAperiIpgOdW1PmpqfVRXB6iqcTan7g9cldV+ysu9Clwik0DBg5Qxphy4DieMrLfWbjr6fWvtvcC9xpivGmOqrLX/WYh+nsRGnKKcw20l0/9eQfbOG9Kw1u0AACAASURBVKm+tgRVdYPrPDJ7E3haPMe1+/k9z/Dwg1upqSqmvbuVDS+u5qLzf8Sic8r58tfeSmXV8JW8RSZSSYmLuvoK6uorTtn2kkucB8ArV82jZkYZKzoOcvjQDwiHiojmM1JF2SFmNzxDSXF84Fxn5CqH+6hwlc46U+NHb9lz8KDzXF8/uF1Pd2gWO/e9ktmNTtWTbDbHhk0L2LZ7FZdddDue/DU7e2YRCjfS0rSBYlea4uIi/uKtF/GBD12K6yQz6aFQjLIyjxbMi4yjggYpY8xynL3tjj4G8AzwY+DO/ik9a+0txpingckUpJ7OPw9XBqH/vfXDtJkQq1atOuEaqeTBJEUzB0NUam+CwHzfcdMrLz63n//8D+euv6a6/9/evYfHddX3/n8v2brb0ki+xIEktkdOCCGljqRAgYAN1kC4lQKSDRR6ciiWyml7yulFIvQHPaflwcjQw2lLnyIZCqc97VNHSi8cWgqalOQQrrbkXBwgBE0ScvEltjSSbN2l9ftj7z0aSaORNNJoz+Xzep55ZrRnz8zyWJr9mbXW/q7HqCi/yM/P3cbL92/nk595Z6z8gUi2e/u7blu0bXxsksuXrrr1tq4wcPl1XHrhCpcvXeHmW37Cc889ythVw/S0U7fh1peNc+svOkVRPc8/7yylU1EBA26Qmpgsx5iZecFsYrIcMBQVjsa2XRrcy1PP3c5et5jp9PQs//aNPXz4v0I0ClvdkfivfQ1OnYK777a8/51fwGLZsrWS7TtK2L6jjOptW2LLF1VXO71b1dvKqaou19xGkRT43SN1AriXuUDyCuAQUA/UAceNMX3MLSKcUWNGbmFNmF+1fKE697on/S1KzfAPhimtjS/EOcLWBfOjBi5d4RMf7WFqqoBNBTMAVGx5gTe89vv8j09/QCFKcl5JaREvvr4o6VCit3zP5UtX2VpRws4FJy/+1V/Bn/wJjI5OMjY2BTjzCwsKZuYNM45PbqG48CoFZq6Y1sRkOYWbx2J/fwBjY6WUljrhzPO1r8GXvwwf+cgYMzOzAHy/t5Znnt/PG+/4HMZ9znMvvISRqzu4cfc/x167uLSSbdtK2HlNcSxwvemtv8B1N6xt4XKRXOZ3kMJae3jhNren6j1AM04QqcOZdN62sa0DYNsy93cCzcaYBmttol6nw0DEXRQ5Iw3/YJiSV81l1PG+EXYc3hH7eXp6hk9+4mt867tvZHZ2E/tf+lXKSoYpLNrEH/7xL1NRmaAQj0geil++J5GSEqcW1uxsIf+7q9lZK/HSFS5dusLg5Vdy+ZLTw/Wi637I+XPTzM5lJqoqn5s3LwtgeKSEa6+dP9fr+eedocPoQPwcrnKKCkdjIQrg4uV9XBrcw017HoxtO/XQq7hw+SZCr/7z2Lafn/sFxiegvX3+a2zeDFdHXuBvv/QdqqrLqaouo2qb82+vjvu5uHhta3WKZDq/g9TlRButtWeAM0CbMWYvzjyk0+72jRJccJ2QtbbFGNOAs35eTfx9xphWnF40X9bOW2ipyeYjp0fY0Tx3+tLYmSts/czcP/tvv/Qdfnz2Wa675lF+1H+I7535AK+p/d98pPU13HTzrg1rv0iuKCgw7Lymgp3XJF+CaWJiisHLV7nsLkrtLFT9CgYuO7d/7UNjbKkom/eYUMipNu/NwQInSC0MYUtvuzJv23e/t4UnnpgfpH77t+H734d7/n6Q7z34MwAiz7yCickZXlrzL7H9roxWU1a2ieuvn54LW9VO2Jq7Xca+m67RPC7JWmkPUsaYiiSlC4aWuR9r7ZM4Q4BpZYypBbpw5jTFp41G43yNi+KcpVeX4Oy7OqDLGNOP02sWxQlPh4HQEj1VGcHOWCZHZigodeZG2MlZCqZmKQw43yIvnBviH08609iu23WWii0XOX/pRt7xrr286W2qkSOSTsXFhex6UWBVZ8L+9m871088XkrDnS9j4PJVrrvhW1y8OMXk+Nx+O7c9gbXz50RNTpVTXDg/XA0MbI6Vb/CcP+/0eg3G9XpdGtzD5NT83ukf97+B8YmtvJYv89wzTvX6n0QOMjOzmZfd+E3A6U27+09+l5JSuP32xf+eH599nqefuuSEr6oyAtVlVFWVaz1FyRgb8ZsYNcYM4ky27sE5M+8p975jOOHlTRvQjqSstX3GmEQhCXAWHHb3W3S/uy3khrF6nF6sLmttRtX2TjTZ/MqjVyh66dwwxPhjV6mImy/1t3/9HaanZ2M/V2y5SF2d4b985L1pb6+IpO7Gl+zi9z725nnbpqZmiA5cZWDgKoOXrzJw+SoDl1/F4MAoA5evsO+mHl54YZKrIwVMT89SXLKZu5oMxQumQA4NOVXmBxYMH8YXQgWnfERR0ei8bQND11Fg5sYsKwNltLYVMD3tTJL3vOc9zrqdd7zicf65qxeAC5f2UbBpmh1VT1FWXkRVVTmVVWXzAtbcz+UEAmUEqsq0rqKk1YYEKaAKp4emEcAYE8UNVUDYGHMSOOp30c1kdZ5WUgPKXYi4b10blWaLCnGeGaHSnWj+ZP8LnPyHUZ49/zZu3P0dyssGKSgw/N4fvlnfBkWyUGHhJnZcU8GOZYYUZ2ctV0bGGR4e47rrFweQs2dhchL++q8mY9v2XHeawk0Ti/YtKZofriYmywlsPRf7OVBVxoVH4ZZb5j/uzBnYs2d+r9cTT99BUeEoO6qeYvTqJKNXJ+n+2lsoKRrhtlu+Gtvv6eduo7BwjBft/AkAv3jb9fzBJ45QUeGEs4V+9vgFTIEhUFVGZaCUzZt19qKs3EYcDQdwep6iwBtxzsqrwhn2ip871GCM+RRwb1yPlaTZyA9GKPm1uXlOY71XuP7jNwBOb9TI1W1cuHwjFy7dxM013+JDH5pi957tfjVXRDZAQcHcUj5LKSqC3/idN/DBD7+OocFRBi5fZXBwlMGBPQwOXGXg8hVe87rvudsCDFy+wsT4NDcH76cortRDVXU5u3bBzTfPf/4LF+CVr1w4ab6MLeUvzNtvfGLrorD25HO3U1F+MRaktlSU8MpXwk03wb/929x+H/kI3Hgj/OTMN/nZTy8AzhJAVYFZrr12hkBVGYFAGZVVZXO3A6Vu4HK2q06XbESQiuAMcz2FO9fJnUBey1ywCuKEq+M4JQ+iJB4KlHV25UdX2Rqc+7CcfPwq5bc6Q31Hf/MgJSXf4V//71/zxFOvYXv1Bd7/n9/mV1NFJAMVFW1eUS8XOMtLOcHqKtGB/QwOXCVQVc6n/ufife+9F6qr4Qufc4YGrYXA1vNUbrkQ28damJycP3xoLUxOlc7bFqgq4+JFuOOO+a/xN38Db3kLmLG5fXvPvpvqymfYP/212LyuU4++m6qK59m3+77YfudeeAllJUNUVV5ga0UJ7/nAL/ErTXUkYq3l509fpjKg4JWLNiJINVlrh+I3uBPIn8SpIYUxphLnzLwQc5XCEw0Fdlhrv7UBbc5Jic7a++aN92HcZSbstGVziaGg0Pkjv/ZFAVo//lYa33uRr3R+m+tvuHbZs4xERJZSWlZEaVkRL3rx8pPnDx1yrn/3o3dy6YURBgdHiQ6cY3BwiqHBG51Advkqb3jN3zE9OdfDZW0BL9r5Y6orn4ltKy/fysgIXBNX12tqCgYHYccOy88e8cKacUJY4fxgNjB0A2Ul8w5jnP3pnVy748dUbv0mQ9ExjDHs2gXvfS987nNz+33847B//zRf/LOvAM7k+oLNO9mxE7ZvL6IyUEplpdPTVen2eFVUzr9dVKSpFJks7f87C0NUkn3uJXmwOuz+vFxdJ1mhUluKqZj7FZh6boLS3U7v1AMPOF3eL3oRBPft5I+PvztW3E9EZKO85JZreQmL1/2MNzExRXRwlOjAKNGod72d6GAd0cGrvPTWHVy5ArNxH2EzM/CpT8Ev3DrFT3qdye/T00UUFY7NKwsxPVOMtZvm9XBNzxQyM1s4b1tpWRkXLjj1tTxTU/DJT8KHfn3uhWdn4d//4/1cf+3D3LLP6eGyFnofeze7tj/OdbvOxva9cGkfW8ous337VSoDpbR94m3cfEvi92JsdJKBy1epCJRSXl6sdRg3UEbG3CWCVb2vjcoBC8/aG3lohP6/eD7289TT45TWlDI1Be97n/Nt7Xd+Bz76UaisRN3RIpKRiosLV7ymoqekBO6+G8ZGYWbiDQwOXmUoOsbrXt/jhLJogOjAVWauzLD/5q9SXjYQe+zszGaqKp6hvDRum3XKy++Yq2XMZbdS4pbyuUn5U9MlWArm9XpNzxRxaXAvlVvOx20r5MyPf4W91/2A8rJvMzY6ycxsATt2wB/8AbS2uvtNwx/9EVz/4kv8yz/8PQCmwLCpcCc7d1iqtxVR6fZwefPevJ6urRUlVFaWUhEoo7S0UGc2pigjg9RCbrC6b9kdZVXGI+MU3jC3ENjU0+NUBEsoLIR//Vf42Mfgf/0vZ7X6ypV/PomIZI3SsiLe0bj0Kl8TE1MMRceccDU4ylDUu0SIDhqig3sZGhpl342lPPIIbI87F6egAH7jNyAYHOEnbjnpqalSwFJUODccOTnlFFUtjAtXXk2ueYFrqoxLl+a37/Jlp2ftrl+bO9NwYqKY/3jgA+y74Tvs2/0957HTRZz58Tu4btcjXLvjcQBmZzdx4dI+KrZeoLJihIqKEr5yT/OSay6efz7KUHSMikApFRWlKivhyoogJekx1j9G0Z75Qar0Hc7I6f79ztktTz8NL37xUs8gIpLbiosL2XlNYUrzQ3fudNZXfCpSSFnR7QxFnTD2i/v/nmh0jCsjhYyPTVG4eZxbanqoqnwu9lg7u4mt5RcpLZmrCuSFq/iw5gWr4qKxRfvND2ulXI7uZkf1z2LbxifLefjxt/OSvd+ivLSX0dFJLlzYxE03wZ/+KXz4w85+AwPwiU/AluKnePiUU196dnYTo+M7uGbnBNu2F1BRUer2eJWwtSLuuqKUrZUl7v0lbN1ayqbNuTW6kbYgZYx5l7X2H9P1/H6/Xi4Yj4xTdPtcV9Pk0+OUBEvm7bN790a3SkQkt+wJ7uDXP3wg4X0TE1MMR8cYio4xNLTPuY6OMTQ0ynD0IYaiMDx0HePjU9zyss1873uwd+/c40tL4YMfhG3Vgzztbpuadj7HCxP0esWHq6kFgauispRLl2BsjHlFWJ9/Hv7yL+Hwu+aKcI2OV/Jg7/u5paaHG648HNvW99ivsG/3d9m1/Qn3dUt47sKt7Kh6ki3ll7np5l189i/fz9AQVFXBpgWdXz/4bj9jo5NsrXCC2NaKkozv/Upnj9Q2Y8xJa+2RNL4GAMaYLwCDgILUKoz1j1F+w9xfy5ODl/jkf73C7uD1/Lff28aevVqMWEQknYqLC9lxTeGKykcA/NIvzf85GIQvfQl++L0Krn9xLSPD4wwPjVJX/38Yjuv1KikZ5pZ935xXDBWgYsv5WEX6SjdIQeJeL8PcOoxeCCssnFt3aHKyjCujO5idnUtHY+MBHn/yIEWFo2wpv0xFZQmPPgp1dfDlL8Nddzn7PfGEMyd3S+GTjAw+5DzfVAmXBvdQVfkc5WVX2Lq1hMjP5tcRywRpC1LW2hPGmH3GmCeARmvtw+v9GsaYPThLzISttXev9/Pnusmh6dgaewCPFJzjwe+9hn/5+nX8xV/O8vY3/5ATX75ZJQ9ERDLcK15VwyteVZPwvomJKYaHxhkZHmMo+hKGh8Zil6HoWYaHyxke2s2Lr6vi1a+GRx6BG26Ye/yuXU7x0qELl/FOw5+Z3czmTePz5nBNTicaUiyZt21rRWlsEv62uHPwn30Wvv51aHidiQWTK6PbeOTxt7H/5q9SWvxThqJjGXn2eFrnSFlr29yFfM8YY+5hnepAGWP2Ax8D3g20WGu/uNbnzAfxdaQKbAFfu/nfY/dND0zxqLlA/a3dDAxdx/kXXsLo8I8JBJaehCkiIpmvuLiQHTsL2bFz6/I7A7+wYD36m292amPdH97Lc88EGBoaY2R4jDcOfY2R4UmGhioYGRqjcst5bnvpP1NRfjH22MLNE+zc9gSlxc5cr4qKUgbckx3jg5QXrmamh/BW6Jnr9ZoLZpl49vhG1JHqNMaEgQ7gPjdYfQE4A5xeyfp6xpgKnPIHtcAR9zoM7HOLe8oqbbfbKdk9t1jxyLNXmGAGY2Bb4Bm2BZ7h9aGXak09EREB4GDDS5PePzk5zZVhZ43G4eFxRobG3NA1wPDQXkaGd/HSW1/E694Ad94JZWVzj92/H/7szyxfPTlXub5g0zQVW87Pq+uVl0EKwFobAULGmFrgbuAzgIVY1fIIzpp8Ufe6GgjEXQfdpzJAN/BRa63KIaxSfB2pgfAAz31jMHbf+RcWr8n8ujfcvGibiIhIIkVFm6nevoXq7VuW3XdhSZ19+5xSO69+5Z3OPK/hMUaGxhgefoyR4WsYGa5gZHicpy4mWHXaZxva3WCt7QOa3AKbIfdS514S8YLVCebW3Vu2UrosbzwyTtHuuTP0vvuTq3zvzK+ya8fj7Nr+U7ZtG6f2dp2yJyIiG6OgwFD/yr1J9zl48Iv85PENatAK+TJu44ahbvcCxKqXg9MLNaDAlF5j/WMU3jG33lXk6Ummpkt4/MmDnLv4Uj5+d6/WdxIREVlGxhwp44KTAtQGGI+MU/Grcz1SW8sivLb+2wxf3cn0dDGve4NW5BEREVlOxgQp2Vjj5yao3lEU+/mFC0MYA5VbLlK+pZjb6vf41zgREZEskXnT3yXtrLVMT9rYz+fPTTM4MFdo7dWv3bfkWksiIiIyR0Eqj3h1pK6vvJ5N2+bOfGhsLOT/nf4gP33qNYxPlvPa17/Ex1aKiIhkDwWpPLRrdheF7mLFMzNw476rFG6a4Onn6ikpLue2ep2tJyIishJZO0fKGPMNa+2b/G5HNvHqSF08eZFLEWd9pE2b4Fff81PM+ANYa9hfewObN2tYT0REZCWyLki5Vc5vx6l0LikY6x+jKG5B4osXnBMljbHsrdm+1MNERERkgYwOUu6ixLU4wakWJzwFkjxEVmAsMkbpgSrn9phzxp5nb80Ov5olIiKSdTIqSBljPoRT5dxbVy92V9ztPuBJ4F0b2LScMv7kOBXXFwNw5MgmHjv7ZspKLnDj7u+yN6ggJSIislIZFaSA32AuQEVwQtMpnAWOP22tjQ3nGWO+ufHNyw1TV2cwhc55BqGQ5Zmnnmf4ynaMMVy/p9rn1omIiGSPjApS1tp6Y0wH8IS19rPx9xljLi/YvWXjWpY7ZidmsUVzJ2sePTzKjh88yvgdmxkfL6G4OPMWhBQREclUGRWkAKy1LcaYQ8aYv7LWfjjJfk9uZLtyxfjT4xTeMLc0zOyTE7z15pdx0/+4ycdWiYiIZKeMrCNlrb0P+Kgx5gvuhHOYP09KUjTWP0bhbidIXb4MV/vHKQ2WLvMoERERSSQjgxQ4ixhba38DaDLGHALsco+R5B588EF+612/RZEbpD72sU3Utb2YN/3FtXR1+dw4ERGRLJSxQcpjrf2Me3Obrw3JEbtmd8V6pN7yllnesXuAyipDcbHPDRMREclCGR+kIDbUd8gd6tsPYIz5fZ+blXXuuOMOPvjmD8aC1DveYfn9qif59rctv/zLPjdOREQkC2VFkIJ5Q30hY8w9QLvfbcpGEy9Msaly7hwDe3WGzVsy7pwDERGRrJA1QcrjDvXd43c7stXMrDPV7Px5eC4yzVSVpp6JiIikKuuCFIC1thunurmsgp20bN7pTIb6whcKuP2OUhpOv5z9L/se33/wZz63TkREJPtk85hOm98NyDazU7MUVGwC4C1vsYyPXeYbX7/E7PQgWyqu97l1IiIi2Sdrg5S19l6/25B1ZqGg1AlSr3iFZbN5gYtPfQ2AvcFDfrZMREQkK6V9aM8Y8w1jzDFjzBtSeOz+uIKcskZ2xmJK5/7LJyemAdi5q4LyLap/ICIislobMUfqdqAV6DHGzBhjTq0iWG0DjruP+VB6m5n77KxlthzOnYNHHoErV5wgtSe43eeWiYiIZKeNCFIDwH3AGZxlXupw5jctG6ystfdZaw8DDUCVMebYBrQ3d83AxFbLv/xLAXfeWch/+vVf5NGfvom9NTv8bpmIiEhW2og5UlGg1Vr7EIAx5jacYPRG4BBOsKoFWo0xAH1AGOix1v4HODWkgM8YY76xAe3NWXbWMrHV8qYDs2zdavn3f32Bc8+cZ29QE81FRERSsRFB6hgQ8X6w1p7B6Z36DIBbXLMBOO1e17mX+GAVAQIb0NbcNgPTW+DGvbB3r6W86Gd8tfth9tTc5nfLREREslLag1Sys+vcobofusN33rZEPVZ17t21aWxqzrOzlpnyuZ8nJ6bYvKmA666r9q9RIiIiWczvgpwN1trPxm+w1p6x1n7GWhuy1hbgBKoIUOsND0pqnoo8xcDEZr75TcOzz8LE+DTXXVPFps1+/xqIiIhkJ7+PoGa5Hay1YaAeOGGMqUh/k3KXKdjEjx4v5667NvOKVxTyF194Bbuv3+Z3s0RERLKW3wU5I8aYd1pr/ynZTtbaqDGmBTgBHNmYpuWenTuv4WUvm+Fzn5tmdNTQ+4Mn2HmNsqmIiEiq/A5SnwZOGWOarbVfSrajtbbPGBPcoHblpNkZy96g4bUHLGAZvvhzSrdc63ezREREspavQ3vW2j7gszjDdj80xvz6MsN3yw4FytJmZy3FRYWxnycnpyndWuRji0RERLKb3z1SWGvb3DIHf4Bzdl6nMaYPpxxCL+5Ec5whvVN+tTMXzM7OcvFSIRdegLIyGB4poERBSkREJGW+BymIhakeoAPYy1yRTo8Beq21H/ajfblidtbyN39Tyt/+rdMrdevNt/CBuzb53CoREZHslRFBCmJn59UYYw4BTUAQqMbpkeqx1p7ws325wFpoarLcdpsz2fzB+yKUlN7id7NERESyVsYEKY+19j6ctflknW2ZLaS+voD6emey+avu3cb+uhv8bpaIiEjW8ruOlPhoy2whxcWFy+8oIiIiCWVcj5SkkYWf/hSmppzJ5tNW86NERETWQkEqj1jgd393E319TkfkJ2sCvM7fJomIiGQ1Bak883u/N8u5c7OMDltu/fq4380RERHJagpSeeb1r7cATF+eYvDUtM+tERERyW6abJ5PCuYKw9uxGQrK9N8vIiKyFuqRyiPjRdM89JChuNgyfWmap6ee51Zu9btZIiIiWSsrgpQxZg/QhlOkE5winb3W2i/61aZsNFk0y1vfuglrDVDFm155jqN+N0pERCSLZUWQAvrd64i19kZjTCXQYIz5Jk6gutvHtmWPAsPnPz/D6Cj8/Olhzj064XeLREREslrGT5JxQ5NxL9UA1toha+291to3AhFjzEk/25gtTAG8852WX/1Vy9vfOsDE1Ut+N0lERCSrZXyPlLV2yBhTBzQA3QnuPwFoHb6ViJtsPjk5DTMzPjZGREQk+2V8kAKw1p4BzvjdjmxnjaGvz1BWZnn2GcOe3VpnT0REZC2yIkjJ+pia3sTb3ub8l5eX1/Bf3nba5xaJiIhkNwWpPFJUBJ///DRjY/Djs89QbLVgsYiIyFooSOWRwkJ417ucyub/9s/nOfeQgpSIiMhaZPxZex63lpSsgTFxk82npikpUZASkdwWCoWoqqrCGENnZ6ffzZEclDVBCqgxxpw2xhxTqErN7KzhzBnD44/DCy8YSssUpEQkt/X09NDc3AxAfX29z62RXJQ1Qcpae5+1th4YBPqNMf/ud5uyzchIIW9962Ze//pC/u4fbqZYQUpE8kA4HAagtrbW55ZILsq6OVLW2uPGmG3A7/vdlmyzdessn//8NKOj8EjvU5SWK0iJSG6LRqP09fXR0NDgd1MkR2VUj5Qx5tgKdz2GU+lcVqGk1Jls/v73W/bufo7SLcV+N0lEJK283qhQKORzSyRXZVqPVAOw7Lp51tqoMSayAe3JKQXxk80npinZqh4pEckdfX19tLW1ARAIBDhy5Ag9PT0ANDY2+tk0yWGZFqSCxpjd1tqnV7CvgtQqTU4aHn7YUFpqGR4poLRCPVIikhu6u7s5evQoXV1dsWG8lpYWOjs7CQQCBINBn1souSrTglQVziLE/Tjr6oWttf+xxL5245qVG85f2MSb3+z8l9fv38ev3VXkc4tEZD088uZHGHtyzO9mpKx0bykv//rLU358OBymqamJjo6OeXOhvCCVbH5UZ2cnHR0d9PX1Ac6E9CNHjtDa2ppyeyS/ZFqQAmdNvVrgo0CbW/uoFzgJ3GetfcjHtmW1a3bCH/3RNKOjhtPhZynfeqvfTRKRdTD25Bhjj2dvkFqrlpYWAoFArMyBJxJxBi6SzY9qbm6mvr6euro6Ojo6Fj2HyHIyLUj1uSUOMMYcAt4INAL17sW6wSrs/iyrUFEB7363BSyv+/vt/ML+6/1ukojImoTDYSKRSMIA5M2PWu6MPW9Cus7sk1Rk1Fl7wFHvhls3qs1aW4Mz5NcEfBEYAkJAwJ8mZrG4yeY65VFEcoEXlhL1OoXD4RXNj+rp6dE8KklZRvVIWWvPLLF9CLjXvbQYY4LAqY1sWy64ehUefthQVmaxM5mWoUUkVaV7S/1uwpqspf3RaBRgUQiKRqNEIpHY2XqRSIRoNJqwKGc4HNZZfZKyjApSK2WtjRhjev1uR7Z54gkTm2z+2WAFB3xuj4isj7VM1M52NTU1gFPuIJ63rp7XU9Xe3h4rjRDPm2SuOlOSqmzullj8FyFJ7dlj+fM/n+bYsWlesnXc7+aIiKxZc3MzgUAgNs8JnBDlDflVV1cDMDAwkHDobrn5Ud6EdZGlZGWPFCw9DJgtjDEBoB3w/rKrgdNAu7U2LX+5gQA0NlpmR2d44Rsz6XgJEZENFQgE6O3tpaWlJRaevEKcx48f59ixY5w8eZL29vaEjz916lTS+VHxzyuSSNqDlDHmG0Af0JOkJtRSgEYLZwAAIABJREFUj90PRK21T6WjbX5xQ1QX0OKFprht/caYOmttX7pef3ZsloKybO6MFBGZEwwGE4ad1tbWZetBhcNh6usTnwTe2dmpIT9Z1kb0SN2Os/RLq1u6oA+nfMFKgtU24GPGmL1Ah7X2i2lt6ca5m7meKCC27E0TMAicAOrW+0WHh+GRR2D26hQPTfXzcvJ3XoWIiDcBPVFY6uzspKWlhcHBQR9aJtlkI4LUAM6QVTVOoc0697JssLLW3gfcZ4ypBJqNMcestcuuxZcFLuO8H/O4YSqK8z6tu8ceM9x5ZyFQya+8OpqOlxARyXjRaJSmpiZOnz4NQEdHR6xHa2BgIBawGhoaFk1iF1loI4JUFGj1KpIbY27D6aF6I3AIJ1TVkiRYueUPPuMOE2YkY0wtzpynapzwWI3Ti9a5cF9r7XHg+BJPFcB5D9bdTTdZfvM3p3nm5yOc75tOx0uIiGS8QCCgeU+ybjYiSB0jboFhd5L4GeAzAMaYe3CC1Wn3OlGPVQQfCnAaY9qBU9ba7mX2awY6gJC1NuxuC+D0pjVZa1c0yG6M8Qbzj62h2Uvatg2amixP/GSYvzurkpwiIiJrlfYgZa29d6n7jDHHgB9aaw/HbUvUY+XNF0rLkNeCNgXd12/DmceUtMyC2xPVgTNxPHb+rTtMdwgYNMa0W2uXe55GnB6ttuWC21pNTExRXJi1J2yKiIhkDL9P3Wqw1n42foO19oy19jPW2pC1tgAnUEWA2nQuWGyMaTTGDOKEogBxvWjLOAGwxBBeFGeYstUNaEu9dgNOiAq5w35pMTwMjz4KTz9tKFKQEhERWTO/g9Sy40tuL089cMIYU5Guhlhru621VW6AO44ztyspd/iuluRzmryB+ITrD7g9UW1AXdywYFpWznzsMcOb3lTIR//wBoqLFaRERETWyu8gFTHGvHO5ndyenRbc3p8M4gWegST7eD1bRxbe4c6Jut0Nb/HBLS1V22+6yfK5z03z3sPPUVJcmI6XEBERySt+d0t8GjhljGm21n4p2Y7W2r5kw2M+8SaRJxsG9ALSvPldxpgO4DBwj3sbnDP90japfts2OHLEEv76C5z7oYKUiIjIWvkapNxw9FmcYbsWnPlJXdba4SUekmmnmnnBLtkwYKy3yhgTcCehtwLN7ubmBI9J2zwpgMnJaUpKFaRERETWyu8eKay1bW6Zgz/AOTuv0xjTh1MOoRd3ojnO0Ngpv9q5hEVFNVewf3SZOlJp4002v3CxgIe+8wD/eDDx2lN33XUXd91118Y2TkREJAv5HqQgFqZ6cHqk9jJXpNNjgF5r7Yf9aF8SKxmGi++t8rVErjfZ/OW37mVL9Od899kHEu538ODBjW2YiIhIlsqIIAWxs/Nq3NpLTTjDZtU4PVI91tpMm2iedW680fLhD0/zo4d/zuUfb+dAzYGE++3Zs2djGyYiIpKlMiZIebz19fxuxzqK74XydYG77dvhPe+x/PXgRd5xx3/ibR9a93WRRURE8orf5Q+y3WqDUbIyCRtmYmKKki2abC4iIrJWClJr4wWjFc19WlArasONjMBjj0E0WsDv/P5va1VzEcl5oVCIqqoqjDF0di5agEJkzRSk1saraJ6svpV330qXnEmbs2cNoVAhZx65gYnZUb+bIyKSdj09PTQ3O1Vm6uvrfW6N5CIFqbXxyjEkK4Pg3RdOss+GuPFGy5/+6TQvr7hE+P/9G9Gorx1kIiIbIhx2Pn5ra9O+7r3kIQWpNbDWdrs3k/11ejO6e5LssyG2b4f3vtfygYsV7Nm3w+/miIikXTQapa+vj4aGtCxhKqIgtQLblrm/E5IuNHwYiMSFLv9ZvxsgIrIxvN6oUCi0zJ4iqcm48gcZJLjgOiFrbYsbojqAmvj73KVgAjh1sXw3MgI/+hEMXBjgTneiuYb3RCRX9PX10dbmrPkeCAQ4cuQIPT3OYEBjY6OfTZMcpiDlMsbUAl0sXji40RhjcUodDAB1Cc6+qwO6jDH9QJu7bxNOb1TILTbqu7NnDQ0NhXx8a6n6IkUkp3R3d3P06FG6urpiw3gtLS10dnYSCAQIBjNtzXvJFQpSLncB5UQhCXAWHHb3W3S/uy3khrF6nF6sLmttSzrbvFo33mhpaZnmlr+p5I8fVk+USK756r19/N9/emhNz7Hr2gr+5DNL9950fv5bnPr+k2t6jf21N/Cbv7t+c5bC4TBNTU10dHTMmwvlBalU5kfV1NQsej6RRBSk4iSr87SSGlDW2j7mSiJknO3b4X3vszzfNe13U0QkDYaHxnj25+mt+3v50pU1v8b1N6x2vffkWlpaCAQCsTIHnkjEqTqz2vlRLS0tRCIRenp6FKRkWRrgyTOzozNsKtvkdzNERNZFOBwmEolw+PDhRfd586NWE4a8yengzLkSWY6CVB4ZGYEzj03zf+xZv5siIrIuvLCUqNcpHA6van5UNBqlra2Njo4OgsEgp0+fXte2Sm7S0F4eOXvW8PZ3bOG9d4z73RQRSYOKylKuW+Ow2a5rK5Lev237ljW/xrbtW9b0+HjemccLw1I0GiUSicTO1otEIkSj0aRFOY8ePcqJEycAp3hnd3c30WhUy2lJUgpSeWTfPst7jgzz2L8OxD4YVP5AJHf88rtr+eV3p7d6d/NvvZ7m33p9Wl9jNWpqnKozC8OOt66e11PV3t4eK42QSGdnJ8FgMBa0br/9drq7uzl9+rTmSUlSGtrLIzt2wFvfMkKB1WRzEckNzc3NBAKBeXObOjs7Y0N+1dVO79nAwMCSQ3yRSISOjg7a29tj27xApXlSshz1SOWZyYlpdl9/PdGH1BMlItkvEAjQ29tLS0tLLDx5hTiPHz/OsWPHOHny5LyQtFBTk1MzOdE8q1OnTi3aJhJPQSqPjIzAz/qhuFD/7SKSO4LBYCxExWttbaW1tTXpY48fP05LS8ui0gkAVVVV6pGSZWloL4+cPWv4o/++i+JiBSkRke7ubk6dOpUwRIFTNsGrRSWyFAWpPLJvn+V977lASXGh300REfFNX18fNTU1NDU1EQ6H582v8oRCIbq7nbXma2pqYpPXRRZS10Qe2bEDam+7zLlvK0iJSP6qra2lv78/6T6JhgpFElGPVJ6ZnJimJ/zvBAIB1UYRERFZIwWpPHLlCpw7X8DM7KTfTREREckJClJ55NFHDff+03Xc9cEPEI1GVYxTRERkjRSk8si+fZZX3v4spVuK/G6KiIhITlCQyiM7dsD1L7pMyRZNNhcREVkPClJ5ZnJymtKtxX43Q0REJCcoSOWRK1dgcNBQWqGhPRERkfWgIJVHHn3U8PDZaymrVI+UiIjIelCQyiP79lleU3KFO9/+etWREhERWQcKUnlkxw6481IpV6cH/G6KiIhITtASMflmepaBIQUpERGR9aAglUeuXIFxu8nvZoiIiOQMDe3lkUcfNfRPlPjdDBGRDRMKhaiqqsIYQ2dnp9/NkRykIJVH9u2z3FA64XczREQ2TE9PD83NzQDU19f73Bp/dHZ20t3d7XczcpaCVB7ZsQMCRTN+N0NEZEOFw2EAamtrfW6JP3p7e/1uQk5TkMozmzbrv1xE8kc0GqWvr4+Ghga/m+KLaDSqIc0001E1j4yMWP5q82m/myEismG83qhQKORzSzZeX18fhw4d8rsZOU9BKo+cPVvAbIHxuxkiImnR19dHKBQiFArR1NREd3c3PT09ADQ2Nm5YG2pqajDGYIyhqqqKvr6+effH3+cFPYDu7m6amppoaWmhpaWFpqYmIpFIwteJRCI0NTXN+/fG7x8Ohzl69CjRaBSAtrY26urqqKuro6WlZcm2e8/jXRL1ZoVCodi/MRKJ0N3dHduWl3OxrLW65PgFuB+we/fcYQ//YrutrKyMXUREckFXV5cNBAK2p6cntq25udkCNhAIbHh7vNeOb4+no6PDNjQ0LNo/EAjYwcHB2Laenp6Ez9HR0WEB29XVNW97IBCwwWBw3raurq6E+y7U3t5ug8Gg7e/vX9Su2traee3q7++3ra2tFrDt7e22vb3dWmstYGtra5O+zlodOHDAAha432bA8dVaqyCVDxcvSNXWvsoeufVTClIiOWZ2dtbOzGTvZXZ2dk3/fi9wdHR0zNve29trAdvY2LjkYzs6Omxtba13cLa1tbWxYLAW/f39S752Y2PjvMDiBaNEoauxsXFeOPH+Tc3NzYv2DQaDi15vJUFqqcDmCQQCi57Xe0x827q6umxvb++Sr7MeMjFIqSBnHrHWsnf3bqKPRv1uioisI2vhwoVpv5uRsmuu2YxZw6yDlpYWAoFArMyBxxvmSjY/qrm5mfr6eurq6ujo6Fj0HKkKBoM0NDTQ3d1NNBqNrW0ajUaJRqMEg8HYvu3t7QAJJ8SHQiFaWlpiz3Hs2DGAhMNz/f39KbXVe66lJuQfPnyYzs5OIpHIvHYDHDlyJHZ7o4ZPM42CVB4ZG4WiIlU2F5HcEQ6HiUQiCQOQNz9quTP2vHlK631mX1tbG+FwmGPHjsXC0rFjx2hra4vtE41GY4EvfrsnEonMK9uQjlIOiQJSvJqamthrL3yfkz0uXyhI5ZGnni6kZH+h380QEVk3XlhK1OsUDocJBALLHux7enpWtN9qNTQ0EAwG6ezsjAWpcDgcuw3EeqoCgcC87enU19cXC2LeZHRJnYJUHrnmmimKSxSkRHKNMc7wWLZay7CeFwQWhiCvp8cbbopEIkSj0YQ9OeFwOG3DUm1tbbS0tMTOfks0JBcMBpc8O2+h+vp6wuHwvDC0EgMDc4vVd3R00NbWRjAYjAXIZK/vDRnma2X45aj8QR7ZUj5NSXH2ftiKSGLGGAoKsvdi1pCkvGEnr2fH4wUXr6eqvb190T5ArDRBuupMNTc3x3qblpqD5fVEHT9+POFzhEKhWGD09u3o6Fh2X5gLmPHbBgYGqK6ujv3sDSkuVbrgnnvuoaGhIW8rwy9HQSqPzM5aSkrVIyUiucMLKvH1mDo7O2NDfl5gGBgYSDh0t9z8qJX2FC3XxkgksmSPTmNjI62trbS1tS0KM21tbYRCoVgIrK2tpb29nc7OzkXB6/jx4/P29fYPBAKx98Ob7B6/T3NzM83NzbS1tS3694ZCIaqrq+nq6pq33Qtm8T1d+cpY5/R4yWHGmPuBAzXBOv74yP/kfZ96nd9NEhFZN5FIJHbmHjhnkjU2NnL8+HFOnjxJMBikvb09YZBqamoiHA4zODiY8LlDoVAshKQqGo2yd+9eent7k87D6uvr49ixY0QikVgAbGtrSxjy4vf1ntP7dyfa1yvOGQwG6ejoSNgO7znje69CoRCtra3z9qurq5tXZNR7zo1YhufgwYM88MADAA9Yaw+m/QVXQEEqD3hBqrLyVfz17/8p7/r/XuV3k0REMkJVVRX19fUJw1JnZyfRaHRRkBD/ZGKQ0oSZPFKxdZLSLUV+N0NEJCN4E9ATzY/q7OykpaVlyZ4qEY+CVB4pKpqmREFKRPJcNBqlqamJ06edRdw7OjpiPVIDAwOxgNXQ0JBwgrpIPAWpPGJnLWVbFaREJL/FT74WWSudtZdHpqehpEJBSkREZL2oRyqPXL1ayPt+7QjPTf00tk1VbUVERFKnHqk8UmUso1Oq+SEiIrJe1COVRyom4YfnH2dTuRYuFhERWQ/qkcojW+u2KkSJiIisIwUpERERkRQpSImIiIikSEFKREREJEUKUiIiIiIpUpASERERSZGClIiIiEiKFKRERCRnhUIhqqqqMMbQ2dnpd3MkBylIiYhIzurp6aG5uRmA+vp6n1vjn87OTrq7u/1uRk5SkBIRkZwWDocBqK2t9bkl/unt7fW7CTlLQUpERHJWNBqlr6+PhoYGv5vim2g0qmHNNFKQEhGRnOX1RoVCIZ9b4o++vj4OHTrkdzNymoKUiIjkhL6+PkKhEKFQiKamJrq7u+np6QGgsbFxQ9tRU1ODMQZjDFVVVfT19c27P/4+L+wBdHd309TUREtLCy0tLTQ1NRGJRBK+TiQSoampad6/OX7/cDjM0aNHiUajALS1tVFXV0ddXR0tLS1Ltt17Hu+SqDcrFArF/o2RSITu7u7Ytrybi2Wt1SXHL8D9gD1w4IAVEclFXV1dNhAI2J6enti25uZmC9hAIOBLm7zXj2+Tp6OjwzY0NCzaPxAI2MHBwdi2np6ehM/R0dFhAdvV1TVveyAQsMFgcN62rq6uhPsu1N7eboPBoO3v71/Urtra2nnt6u/vt62trRaw7e3ttr293VprLWBra2uTvs5aHDhwwAIWuN9mwPHVWqsglQ8XBSkRyWVe2Ojo6Ji3vbe31wK2sbFx1c8ZDAYTBqDV6O/vX/L1Gxsb5wUWLxgles3GxsZ54cT7dzU3Nyds98LXW0mQWiqweQKBwKLn9R4T37auri7b29u75OusVSYGKQ3t5ZEHH3yQQCAw7yIiuWN8HKJRcL4/Zce29dDS0kIgEIiVOfB4Q1yrnR/V0tJCJBKJDQumKhgM0tDQQHd3d2x4DZzJ39FolGAwGNvW3t4OkHBSfCgUoq+vL/Ycx44di7Vzof7+frq6ulbdVu+5lpqUf/jwYbq7uxMOMx45ciR2u7GxMe/OjlSQEhHJEZ/+NFRVwdBQ9mxbq3A4TCQS4fDhw4vu84LQas7Yi5+vFD+vKVVtbW3AXPjxbnvbwQlWXkBpa2tbdOnp6ZkXTtJRziESicwLdgvV1NTMe+14yR6XDzb73QDZOHfccQf333+/380QEVk3XlhK1OsUDocJBAIrPtBHo1Ha2tro7e0lHA5z+vTpNbevoaGBYDBIZ2dnrNcpHA7HbgOx0YFAIDBvezr19fXFglh8b5msnnqkREQka3khYGFY8np5vN6oSCSybA/T0aNHOXHiBOD09nhDcGvV1tYWq+XU2dmZcEguGAyu+LW8Cu2r7TEbGBiI3e7o6Ij1gnlhc6mzA8EZMox/bZmjICUikiM++lEYHITKyuzZtlbekNPCOZ/eKfteT1V7e3vSeaGdnZ0Eg8FYL83tt98OsC69Us3NzbHepo6OjkVzubz2ARw/fjzhc4RCoVjQ8vbt6OhYdl+YC5nx2wYGBqiuro797A01LlW64J577qGhoSHv5j+thIKUiEiOKCmBQACMyZ5ta+WFlPi5O52dnbEhPy8sDAwMLDnEF4lE6OjomDes5gWG9Zgn5bUzEoks2aPT2NhIa2srbW1ti8JMW1sboVAoFgRra2tpb2+ns7NzUfA6fvz4vH29/QOBQOw98Xra4vdpbm6mubmZtra2RT1ToVCI6urqRZPYvWAW39OVj4xd79MnJOMYY+4HDhw4cEBzpEQk50QikdiZe+CcRdbY2Mjx48c5efIkwWCQ9vb2JYNUXV0dwLweGnDmMjU2NqZ0FtxC0WiUvXv30tvbm3TOVl9fH8eOHSMSicTa09bWlnDCfPy+3nN6//ZE+3rFOYPBIB0dHQnb4T1nfO9VKBSitbV13n51dXXzQqb3nOleiufgwYM88MADAA9Yaw+m9cVWSEEqDyhIiYgkdvz48YSlEwCqqqqorq6OzQ8S/2VikNLQnoiI5KXu7m5OnTqVMESBc8ZdsgnYIqAgJSIiecZbC6+pqYlwOJywNlIoFIrNVaqpqUm43pwIqI6UiIjkmdra2mWH69Za1Vzyh3qkRERERFKkICUiIiKSIgUpERERkRQpSImIiIikSEFKREREJEUKUiIiIiIpUpASERERSZGClIiIiEiKFKREREREUqQg5TNjTNAY02WMWfvy4iIiIrKhtESMT4wx7UAt0Ac0AosXexIREZGMpiDlE2ttm3fbGJN46fF1YIy5C9iTrufPd1/5yld46qmn2LNnD3fddZffzckpem/TR+9t+ui9Ta/z5897N3f52Y54xlrrdxvynjFmEDhtrQ2l4bnvBw4AHDhwgPvvv3+9XyKvHTx4kAceeEDvbRrovU0fvbfpo/c2vQKBAENDQwBD1tqA3+0BzZFaN8aYWmNMjzGmN+46bT1NIiIi4j8FqSSMMe3GmMYV7NcM9ALt1to6t2fpENBijOlJdztFRETEHwpSC7hn0TUbY/qBViC4zP61QAfQYq2NTRi31kZxwlSDO7FcREREcoyClMsY0+jOVeoAAkBkhQ89AWCt7Vx4hxumwkCrMSZpIBMREZHsoyDlstZ2W2urrLUha+1xILrcY4wxAeZKGCzFG9pbdohQREREsouC1No0uNcDSfbxeraOpLktIiIissEUpNbGK1eQbBjQ69mqTXNbREREZIMpSK2NN+8p2TBgrLfKHQr0bje7S8P04szJanBLJnQZYxoSPZGIiIhkFlU2X5vqFPaPQmxy+qIJ6un04IMPEggsX79s165d7NqVMUVjM9pDDz0Uuz548KC/jckxem/TR+9t+ui9Ta8rV654N0v9bEc8VTZfgruIcCPQ5k4+T7RPP06v1PH4JV8W7BME+t0f66y1ySamrzt3iZgvAMUb+boiIiJpNGatLfO7EaAeqZxnrf2KMeZ3gJcAm4CxFTzsvHsRERHJJPuALcCTfjfEoyCVfvFjacuWVEgHa+1tfryuiIhIrtNk87VZbTBKViZBREREsoyC1Np4wWhFK1C7lc5FREQkRyhIrY03cTzZ8i/efStdckZERESyhILU2pxyr5OVQfDuCyfZR0RERLKQgtQaWGu73ZvJqpbXudc9SfYRERGRLKQgtbxty9zfCZCkGvlhIBIXukRERCRHKEgtLbjgOiFrbQvO/KeOhfcZY1pxJqK3rHvrRERExHeqbO4yxtQCXThzmhKdhRfFOUuvbuHZd+4ael04oavN3bcJpzeqyVqr+VEiIiI5SEEqjjEmsFSJAm/B4WQlDNwwVo8TxPoUoERERHKbgpSIiIhIijRHSkRERCRFClIiIiIiKVKQEhEREUmRgpSIiIhIihSkRERERFK02e8GyPpzyzC049TEGnCvO6y1nb42TER85xYKPoLz2QBOQeG2ZKVdJDn3M/du98cAzmfuSWvtcf9aldsy6TinIJVjjDHNOFXWQ14dK7cG1n3GmCZrbcjXBmY5Y0wjzgdmEOePNwJ0KaSml3vwr3FXEpAUuAee+4B7rLV1cdsb3e11Sz1Wlub+bt4OHI0Po8aYdmPMILBXIXVljDHtwKnlllTLtOOchvZyiPtB2QG0xBcDdf+IDwEN7i+qpMB971pwqtVX4Rx4IkCHMabXK9oq68sYE2Tum6ekwP1s6MUJUYnCaDDJeqGyBPd9DVlrmxaGJWttG3APTkiVJRhjgsaYZmNMP9DKMsuyZeJxTkEqt5wASNQ74v6ShYFW98Akq+B+A6q11oastRFw3lP3oNQJeN/2Zf11+d2AHNAFRBeGqLjlrQIkXhpLkmsn+e9nO1Crz9zFjDGNbo9dB87vXmSFD82445yCVI5wPxBrgb4ku/W4143pb1HO8XqjlroPnA/M5g1qT15w308d4NfAHXoK4gT+eeIOPBH3WlanHqhZ6k7vSxfOZ7PEsdZ2W2ur3C+nx3HWqE0qU49zClK5w+uWH0iyj/dHfSTNbckp7pDH6bgPxXnc7d4ftuagrRP3Q7OFuaAqqfEmQZ9MdKd7IKvRPJ6UNS81rB/XK5LswC8rl5HHOQWp3OEdwJN1j3oflPp2tDq1OOPuvUn20TfP9XcCOEryD01Jwp1P4i24roP5+juN8/4+6b7XC7UD4aW+hMmqZeRxTkEqd3jffJJ9q4wdkDQxelW2ude1SSbkag7EOnLPJIvo4L9m3u9rFGLzUnrckyN6dfLJmrW51wGg1xjT4d3hvre1LD0lQFYvI49zKn+QO1Z7RlM1KxiTFsCZDNmIM1l3qXkk6sJfJ+6H393xp+hLym53rwfc+WY18aeGe6EKOKShvdWz1vYZY1pwPiPAGeZrYK7HpE7v67rKyOOceqRyx0qSd/wvlHqkVshaG3HnkCQ8sMcPnzD3gSqpa2fum76sTfzfecg9JT9eEzrjdE3cs8fqmPt8DeL0BPYpRK27jDzOKUiJrJ03mTecpMdKVsAbOtX7uG68b/BBEoT8uLP2at3hVFklrwfV/TF+7k6rMaZfpQ9yn4JUfolP5/qmtA7cA38jzgdok8/NyQXqjVpfsfkiScKpd7r43UvcL0twQ9KTOJ+te621NUD8sjBBnLlTGgHYOBt+nFOQyh2r/YXRmVBrFFfMMIozbKJwugbu5Nxjeh/XlfdeZtRZTjmkB+ekiNjfvzt8WsPcfMkAKuGxXjLyOKcglTu8X5gVffPRwWpdePNK6nR689q488xql1tjS1ZtVQcS9ZysnDt5P4hTomMed15lHeD9Ph/eyLblsIw8zilI5Q7v20+y8XjvPh3018g9zTmIQtR6OYGGRtMhWe0zWZsQzpm8S56pa61twvlsDiikrouMPM6p/EHuOOVeJzs91LtPE3nXwF1yowGt6r4u3Hkm3lySRLt4v7eN7sKmAB3ushKSnPe3nuzA4x3go/p9XpUgKztYn8TpbdV7u3YZeZxTkMoR1tpu9yCUbJ6Dd/p+T5J9JAn3zKYWlqgPY4zpVf2j1XF79KqWut+dO9UKdLvf8GWFrLURY0wECBpjapfoPfHWirtnA5uWCyLMFTxdjurLrYNMPc5paC+3dMLcKeQJHMaZGKl5KClw5/G0s3SIWuk3VJGN5J0FudTnQsOC/WRljuEM2S1XNuKIu6+szLZl7s+445x6pHKItbbF/eXqYMGK5O5wVADNQ0mJG5Luw/lAbEgwBFWN897qm2f6aI5JCtxv8WGg3RjTHT+nz/1cCAItGnpaHbeqeRtwwhiTcDkjdy6lvryuTHDBdUKZeJwz1tqNfD1Js7hT8oM43zCjOL9Uh4EmFTpcPfc97WVl6+k16UNzfRljenB6TSJunR5JgTGmC+d9PIZbssP9WZ8La+D2SLXjfInqwTmz7Hac+nLtbuVzWcAtZRWEAAADLUlEQVTt4e/C+RKa6EtSFOe9XDQCkGnHOQWpHOX+ktbj/IL26YMydXFzdFaiTgvtrp37jbOLpT9g23SAWj33ANSAc/A6rd/V9RP3mQtO6Ndn7jKMMYGlekK9sxyT9ZRmynFOQUpEREQkRZpsLiIiIpIiBSkRERGRFClIiYiIiKRIQUpEREQkRQpSIiIiIilSkBIRERFJkYKUiIiISIoUpERERERSpCAlIiIikiIFKREREZEUKUiJiIiIpEhBSkRERCRFClIiIiIiKVKQEhEREUmRgpSIiIhIihSkRERERFKkICUiIiKSIgUpERERkRQpSImIiIikSEFKREREJEUKUiIiIiIpUpASERERSZGClIiIiEiKFKREREREUqQgJSIiIpIiBSkRERGRFClIiYisgTGm1RjT714aFtzXbowJ+tU2EUk/BSkRkRQZYzqAdiDoXnqMMY3ufc1A0Fob8bGJIpJmxlrrdxtERLKOMaYWuA84CvThBKla4Hb3OmqtrfOvhSKyERSkRERS4A7jDVhr+xZsDwJdwCFrbdSXxonIhlGQEhFZJ24v1QmgSUN6IvlBc6RERNaB20N1AqcnSiFKJE9s9rsBIiLZzp1g3qI5USL5R0FKRGQN3LPzmqy1Ib/bIiIbT0N7IiIpckNUKFGIMsbUuveLSA5TkBIRSYExphUnRDUtscsJ4PQGNklEfKCz9kREVsmdWN6DUz8qApyKu90AtAERDfeJ5D4FKRGRVTDGBIBenHlRfe5E864Fu/WhOlIieUFDeyIiq3M3bogCsNZ2AyEg6l46UYgSyRvqkRIRERFJkXqkRERERFKkICUiIiKSIgUpERERkRQpSImIiIikSEFKREREJEUKUiIiIiIpUpASERERSZGClIiIiEiKFKREREREUqQgJSIiIpIiBSkRERGRFClIiYiIiKRIQUpEREQkRQpSIiIiIilSkBIRERFJkYKUiIiISIoUpERERERSpCAlIiIikiIFKREREZEUKUiJiIiIpEhBSkRERCRFClIiIiIiKVKQEhEREUmRgpSIiIhIihSkRERERFKkICUiIiKSIgUpERERkRT9/5RDNcwSNnkFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 562.5x562.5 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# now plot the fractional error as a function of step size\n",
    "plt.figure(figsize=(2.5,2.5))\n",
    "plt.xlabel('$z$'); plt.ylabel(r'$d_L, d_A (h^{-1}\\ {\\rm Mpc})$')\n",
    "#plt.xscale('log'); \n",
    "plt.yscale('log')\n",
    "\n",
    "plt.xlim([zmin,zmax]); plt.ylim(10.,1.e5)\n",
    "\n",
    "plt.plot(z, dL, linewidth=1.5, c='m', label='$d_L$')\n",
    "plt.plot(z, dA, linewidth=1.5, ls ='--', c='darkslateblue', label='$d_A$')\n",
    "plt.plot(z, dLv, linewidth=1., c='lavender', label='$d_L\\ \\mathrm{vector}$')\n",
    "plt.plot(z, dA, linewidth=1., ls =':', c='blue', label='$d_A\\ \\mathrm{vector}$')\n",
    "\n",
    "plt.legend(loc='lower right', frameon=False, fontsize=7)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the peculiar behavior of $d_A$ at high $z$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Peculiar velocities\n",
    "\n",
    "For nearby galaxies (low luminosity) or for galaxies near very large structures distances become inaccurate due to *peculiar motions*.\n",
    "\n",
    "If physical distance between two galaxies is $r$, comoving distance between is $x=r/a$ (where a(t) is normalized to 1 at $z=0$). \n",
    "Thus, relative velocity between these two galaxies is \n",
    "\n",
    "$$\n",
    "\\dot{r}=\\dot{a}x + a\\dot{x}.\n",
    "$$\n",
    "\n",
    "If galaxies were simply coasting with the expansion of the universe, we would have $\\dot{x}=0$ by definition and pure Hubble law:\n",
    "\n",
    "$$\\dot{r}=\\dot{a}x=\\frac{\\dot{a}}{a}\\, r=H(a)r;\\ \\ \\ \\ H(a)=\\frac{\\dot{a}}{a}$$\n",
    "\n",
    "However, because distribution of matter is not completely homogeneous, galaxies feel forces of gravity due to local inhomogeneities that perturb their motion away from the pure Hubble flow. Thus, $\\dot{x}\\ne 0$ and\n",
    "\n",
    "$$\n",
    "\\dot{r}=cz=H(a)r + av_{\\rm pec},\\ \\ \\ \\mathrm{where}\\ v_{\\rm pec}=\\mathrm{peculiar\\  velocity}.\n",
    "$$\n",
    "We measure distances from $z$ as $d=cz/H$. Thus, what we are measuring is\n",
    "\n",
    "$$ d = r + \\frac{a}{H}v_{\\rm pec} = r + \\frac{1}{H}\\frac{v_{\\rm pec}}{1+z}.$$\n",
    "\n",
    "So the positions of galaxies in the regions of massive clusters or superclusters (or any time we use redshift as proxy for distance) will be distorted by $v_{\\rm pec}/(1+z)/H$, which is manifested in the \"streaks\" in galaxy distribution pointing towards us. Most prominent ones are strangely called the *\"fingers-of-God\"* (because they all are extended structures \"pointing\" at us) and are due to galaxy motions inside massive galaxy clusters. "
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}