guyer/Euler.ipynb

## Euler.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f2be82be-905d-4af5-a8d8-1fb0d1c04655",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Smallest workable example of hydrostatic equilibrium. Norbert 22.01.2024\n",
    "\"\"\"\n",
    "\n",
    "import fipy\n",
    "import numpy as np\n",
    "import pylab as plt\n",
    "\n",
    "\n",
    "#system of units (relative to CGS): \n",
    "# user specified units\n",
    "unit_length = 1.e8 # cm\n",
    "unit_numberdensity = 1.e9 # cm^-3\n",
    "unit_temperature = 1.e6 # K\n",
    "\n",
    "# constants\n",
    "He_abundance = 0.1 \n",
    "mp = 1.67e-24 # g\n",
    "kB =  1.3806488e-16 # erg K^-1\n",
    "miu0 = 4*np.pi\n",
    "gamma=5./3.\n",
    "R_sun = 6.957e10/unit_length # cm\n",
    "a = 1.0 + 4.0*He_abundance\n",
    "b = 2.0 + 3.0*He_abundance\n",
    "# derived units\n",
    "unit_density = a*mp*unit_numberdensity \n",
    "unit_pressure = b*unit_numberdensity*kB*unit_temperature\n",
    "unit_velocity = (unit_pressure/unit_density)**0.5\n",
    "unit_time = unit_length/unit_velocity\n",
    "    \n",
    "L = 1.e10/unit_length # in cm\n",
    "nx = 1000\n",
    "dx = L/nx\n",
    "mesh = fipy.Grid1D(dx=L / nx, nx=nx)\n",
    "\n",
    "x = mesh.cellCenters[0]\n",
    "x.name='x'\n",
    "\n",
    "g = 2.74e4 * unit_length/unit_velocity**2 # test with constant g \n",
    "\n",
    "T_ini = 1.e6/unit_temperature # in K \n",
    "\n",
    "    \n",
    "dens=fipy.CellVariable(name=\"dens\",mesh=mesh,value=0., hasOld=True)\n",
    "scaleheight = 5.e9/unit_length #in cm\n",
    "scaleheight = T_ini*unit_temperature/(g*unit_velocity**2/unit_length)* kB/mp/a*b /unit_length\n",
    "basevalue = 1.e-15/unit_density # in g cm^{-3}\n",
    "dens.setValue(basevalue*np.exp(-mesh.cellCenters/scaleheight))\n",
    "dens.constrain(basevalue*np.exp(-mesh.faceCenters[0].value[0]/scaleheight), where=mesh.facesLeft)\n",
    "dens.constrain(basevalue*np.exp(-mesh.faceCenters[0].value[-1]/scaleheight), where=mesh.facesRight)\n",
    "        \n",
    "vr=fipy.CellVariable(name=\"vr\",mesh=mesh,value=0.)\n",
    "basevalue_v = 0.e0/unit_velocity # cm/s\n",
    "vr.setValue(basevalue_v*np.exp(vr.mesh.cellCenters/scaleheight))\n",
    "    \n",
    "vel = fipy.CellVariable(name=\"vel\", mesh=mesh, value=0., hasOld=True)\n",
    "vel.setValue(vr)\n",
    "vel0=vr.value[0]\n",
    "vel.constrain(0.0, where=mesh.facesLeft)\n",
    "vel.constrain(0.0, where=mesh.facesRight)\n",
    "\n",
    "# velocity at faces\n",
    "vel_f = fipy.FaceVariable(name=\"vel_f\", mesh=mesh, value=0., rank=1)\n",
    "vel_f.constrain([[0.]], where=mesh.exteriorFaces)\n",
    "\n",
    "Eth = fipy.CellVariable(name=\"Eth\", mesh=mesh, value=dens.value*T_ini/(gamma-1), hasOld=True)\n",
    "Eth.faceGrad.constrain(value=g*dens.faceValue * [[1.]] / (gamma-1), where=mesh.exteriorFaces)\n",
    "\n",
    "EthPrevious = Eth.copy()\n",
    "densPrevious = dens.copy()\n",
    "velPrevious = vel.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d762b947-6157-43c8-9eb2-537e1c33a0d0",
   "metadata": {},
   "source": [
    "We obtain the energy-corrected velocity by starting with the momentum equation\n",
    "\\begin{align}\n",
    "\\frac{\\partial\\rho\\vec{v}}{\\partial t}\n",
    "&= -\\nabla\\cdot\\left(\\rho\\vec{v}\\vec{v}\\right)\n",
    "- \\left(\\gamma - 1\\right)\\nabla E_{th}\n",
    "- \\vec{g} \\rho\n",
    "\\end{align}\n",
    "following the [FiPy's scheme for discretizing linear equations](https://pages.nist.gov/fipy/en/latest/numerical/discret.html#section-linear-equations), and discarding terms implicit in variables other than $\\vec{v}$\n",
    "\\begin{align}\n",
    "\\frac{\\partial\\rho\\vec{v}}{\\partial t}\n",
    "&\\approx -\\nabla\\cdot\\left(\\rho\\vec{v}\\vec{v}\\right)\n",
    "- \\left(\\gamma - 1\\right)\\nabla E_{th}\n",
    "\\\\\n",
    "\\frac{\\rho_P\\left(\\vec{v}_P - \\vec{v}_P^\\text{old}\\right) V_P}{\\Delta t}\n",
    "&=\n",
    "-\\sum_f \\rho_f\\left(\\vec{n}\\cdot\\vec{v}\\right)_f A_f \\left[\n",
    "\\alpha_f \\vec{v}_P + (1 - \\alpha_f) \\vec{v}_A\n",
    "\\right]\n",
    "- V_P \\left(\\gamma - 1\\right) \\left(\\nabla E_{th}\\right)_P\n",
    "\\\\\n",
    "a_P \\vec{v}_P &= \\sum_f a_A \\vec{v}_A - V_P \\left(\\gamma - 1\\right) \\left(\\nabla E_{th}\\right)_P\n",
    "\\\\\n",
    "\\vec{v}_P &= \\sum_f \\frac{a_A}{a_P} \\vec{v}_A - \\frac{V_P}{a_P} \\left(\\gamma - 1\\right) \\left(\\nabla E_{th}\\right)_P\n",
    "\\\\\n",
    "a_P &= \\frac{\\rho_P V_P}{\\Delta t} + \\sum_f\\left[a_A + \\rho_f\\left(\\vec{n}\\cdot\\vec{v}\\right)_f A_f\\right]\n",
    "\\\\\n",
    "a_A &= -\\rho_f\\left(1 - \\alpha_f\\right) \\left(\\vec{n}\\cdot\\vec{v}\\right)_f A_f\n",
    "\\end{align}\n",
    "\n",
    "We obtain the Rhie-Chow correction for collocated velocity and energy by letting\n",
    "\\begin{align}\n",
    "\\vec{v}^\\diamond_P &= \\vec{v}^*_P + \\frac{V_P}{a_P} \\left(\\gamma - 1\\right) \\left(\\nabla E_{th}^*\\right)_P\n",
    "\\end{align}\n",
    "be an intermediate velocity without energy correction and\n",
    "\\begin{align}\n",
    "\\vec{v}_f &= \\left(\\vec{v}^\\diamond\\right)_f - \\left(\\frac{V_P}{a_P}\\right)_f \\left(\\gamma - 1\\right) \\left(\\nabla E_{th}^*\\right)_f\n",
    "\\end{align}\n",
    "be the velocity at a face, in terms of the intermediate velocity, but with energy correction reapplied. Expanding,\n",
    "\\begin{align}\n",
    "\\vec{v}_f &= \\left(\\vec{v}^*\\right)_f + \\left(\\frac{V_P}{a_P} \\left(\\gamma - 1\\right) \\left(\\nabla E_{th}^*\\right)_P\\right)_f - \\left(\\frac{V_P}{a_P}\\right)_f \\left(\\gamma - 1\\right) \\left(\\nabla E_{th}^*\\right)_f\n",
    "\\\\\n",
    "&= \\left(\\vec{v}^*\\right)_f + \\left(\\frac{V_P}{a_P}\\right)_f \\left(\\gamma - 1\\right) \\left[\\left(\\nabla E_{th}^*\\right)_f - \\nabla_f E_{th}^*\\right]\n",
    "\\\\\n",
    "\\vec{v}_f \n",
    "&\\approx \\left(\\vec{v}^*\\right)_f + \\frac{A_f d_f}{a_f} \\left(\\gamma - 1\\right) \\left[\\left(\\nabla E_{th}^*\\right)_f - \\nabla_f E_{th}^*\\right]\n",
    "\\end{align}\n",
    "where\n",
    "\\begin{align}\n",
    "\\left(\\vec{v}^*\\right)_f &\\qquad\\text{is the average value at the face of the current value of the velocity}\n",
    "\\\\\n",
    "\\left(\\frac{V_P}{a_P}\\right)_f &\\approx\\frac{A_f d_f}{a_f}\n",
    "\\\\\n",
    "V_P&\\qquad\\text{is the volume of the cell}\n",
    "\\\\\n",
    "a_P&\\qquad\\text{is the matrix diagonal for the momentum equation}\n",
    "\\\\\n",
    "A_f&\\qquad\\text{is the area of the face}\n",
    "\\\\\n",
    "d_f&\\qquad\\text{is the distance between the centers of the cells adjoining the face}\n",
    "\\\\\n",
    "a_f&\\qquad\\text{is the average value at the face of}\n",
    "\\\\\n",
    "&\\qquad\\text{the matrix diagonal for the momentum equation}\n",
    "\\\\\n",
    "\\left(\\nabla E_{th}^*\\right)_f&\\qquad\\text{is the average value at the face of the gradient in energy}\n",
    "\\\\\n",
    "&\\qquad\\text{as evaluated in the cells adjoining the face}\n",
    "\\\\\n",
    "\\nabla_f E_{th}^*&\\qquad\\text{is the gradient in energy evaluated at the face}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "aedc68c6-ba83-4ab9-8c55-929a8cf382f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "matrixDiagonal = fipy.CellVariable(mesh=mesh, name='$a_f$', value=1e+20, hasOld=True)\n",
    "correctionCoeff = mesh._faceAreas * mesh._cellDistances / matrixDiagonal.faceValue\n",
    "\n",
    "# uncorrected velocity at faces\n",
    "uncorrected_velocity = vel.faceValue + correctionCoeff * (gamma-1) * Eth.grad.faceValue\n",
    "\n",
    "# corrected velocity at faces\n",
    "corrected_velocity = uncorrected_velocity - correctionCoeff * (gamma-1) * Eth.faceGrad\n",
    "corrected_velocity.name = \"vel_f\"\n",
    "corrected_velocity.constrain([[0.]], where=mesh.exteriorFaces)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3495b55-e3ae-46fd-baaa-a664d18a9bfc",
   "metadata": {},
   "source": [
    "Substituting this energy-corrected velocity in the continuity equation\n",
    "\\begin{align}\n",
    "\\frac{\\partial\\rho}{\\partial t} &= -\\nabla\\cdot\\left(\\rho\\vec{v}\\right)\n",
    "\\\\\n",
    "&= \n",
    "-\\nabla\\cdot\\left\\{\n",
    "\\rho \\left[\\vec{v}^*_f\n",
    "+ \\frac{A_f d_f}{a_f} \\left(\\gamma - 1\\right) \\left(\\nabla E_{th}^*\\right)_f\n",
    "\\right]\n",
    "\\right\\}\n",
    "+\\nabla\\cdot\\left[\n",
    "\\rho_f \\frac{A_f d_f}{a_f} \\left(\\gamma - 1\\right) \\nabla E_{th}\n",
    "\\right]\n",
    "\\end{align}\n",
    "which is represented by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e8bda011-5c8c-4b9b-b9f2-4f4b01e37a61",
   "metadata": {},
   "outputs": [],
   "source": [
    "eq_continuity = (fipy.TransientTerm(var=dens) \n",
    "                 == - fipy.VanLeerConvectionTerm(coeff=uncorrected_velocity, var=dens)\n",
    "                 + fipy.DiffusionTerm(coeff=dens.faceValue * correctionCoeff * (gamma-1), var=Eth)\n",
    "                ) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cc5566b-e55d-41b6-ab3a-68addd4e163f",
   "metadata": {},
   "source": [
    "The momentum equation\n",
    "\\begin{align}\n",
    "\\frac{\\partial\\rho\\vec{v}}{\\partial t}\n",
    "&= -\\nabla\\cdot\\left(\\rho\\vec{v}\\vec{v}\\right)\n",
    "-\\nabla\\cdot\\left[\\left(\\gamma-1\\right) E_{th} \\mathsf{I}\\right]\n",
    "-\\rho\\vec{g}\n",
    "\\end{align}\n",
    "which is represented by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9346af53-8e26-46cc-8cd9-c284a9315a10",
   "metadata": {},
   "outputs": [],
   "source": [
    "eq_momentum = (fipy.TransientTerm(coeff=dens, var=vel)\n",
    "               == -fipy.CentralDifferenceConvectionTerm(coeff=[[1]] * dens.faceValue * vel.faceValue, \n",
    "                                                        var=vel)\n",
    "               - (gamma-1) * fipy.CentralDifferenceConvectionTerm(coeff=[[1]], var=Eth)\n",
    "               - fipy.ImplicitSourceTerm(coeff=g, var=dens)\n",
    "              )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "469c1863-1e78-497a-9158-573ba6080256",
   "metadata": {},
   "source": [
    "Substituting this energy-corrected velocity in the energy equation\n",
    "\\begin{align}\n",
    "\\frac{\\partial E_{th}}{\\partial t} \n",
    "&= -\\nabla\\cdot\\left(E_{th} \\vec{v}\\right)\n",
    "-\\left(\\gamma - 1\\right) E_{th} \\nabla\\cdot\\vec{v}\n",
    "\\\\\n",
    "&= \n",
    "-\\nabla\\cdot\\left\\{\n",
    "E_{th} \\left[\\vec{v}^*_f\n",
    "+ \\frac{A_f d_f}{a_f} \\left(\\gamma - 1\\right) \\left(\\nabla E_{th}^*\\right)_f\n",
    "\\right]\n",
    "\\right\\}\n",
    "+\\nabla\\cdot\\left[\n",
    "E_{th,f} \\frac{A_f d_f}{a_f} \\left(\\gamma - 1\\right) \\nabla E_{th}\n",
    "\\right]\n",
    "\\\\\n",
    "&\\qquad -\\left(\\gamma - 1\\right) E_{th} \\nabla\\cdot\\vec{v}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "eb3e2b28-b304-4336-9b7d-e3eb6d349cbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "eq_eth = (fipy.TransientTerm(var=Eth) \n",
    "          == - fipy.VanLeerConvectionTerm(coeff=uncorrected_velocity, var=Eth)\n",
    "          + fipy.DiffusionTerm(coeff=Eth.faceValue * correctionCoeff * (gamma-1), var=Eth)\n",
    "          - fipy.ImplicitSourceTerm(coeff=(gamma-1) * corrected_velocity.divergence, var=Eth)\n",
    "         )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f40fb43e-2a17-4f1d-9225-d5d01d258ee6",
   "metadata": {},
   "outputs": [],
   "source": [
    "relaxation = 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "046bbd01-da15-4a86-ab2c-6ef52eabbd38",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbMAAAGxCAYAAADoJViMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABahklEQVR4nO3deVwU9f8H8NeywHIvl9zI4Y2IB6aiGSoKeaWZR6moqZX19f5qef3UrL6klZkllqaYpWXmkaWZ5InhkQd5QJ4gHiCCct+7n98fftmvK4uCsiwDr+fjsY/az35m5j0D7ouZ+cyMTAghQEREJGFGhi6AiIjoaTHMiIhI8hhmREQkeQwzIiKSPIYZERFJHsOMiIgkj2FGRESSxzAjIiLJY5gREZHkMcyo1srJycHbb7+N0NBQNGjQADKZDAsXLtTZt1u3bpDJZJDJZDAyMoK1tTUaN26MIUOG4KeffoJara7UMrt16wZ/f/9qXAvd1q1bB5lMhqSkJL0vqzo96mfwOJGRkVi3bl259qSkJMhkMp2fEVUWw4xqrYyMDKxatQpFRUUYOHDgY/v7+vriyJEjiI2Nxfbt2zFr1iwUFBRgyJAh6NatG7KysvRfNFWoojBzdXXFkSNH0Ldv35oviuoMY0MXQFQRLy8v3Lt3DzKZDOnp6fj6668f2d/c3BydOnXSahs/fjyioqIwduxYvP7669i0aZM+S652BQUFMDc3N3QZeqVQKMr93IiqintmVGuVHTZ8Wq+++ir69OmDzZs349q1a5Wa5q+//kLXrl1hYWEBX19ffPjhh5pDlbm5ubC1tcUbb7xRbrqkpCTI5XJ89NFHmrajR4+iS5cuMDMzg5ubG2bPno2SkpJy03p7e6Nfv37YunUr2rZtCzMzM7z77rsAgHPnzmHAgAGws7ODmZkZ2rRpg2+++aZS67J582Z07NgRSqVSsz5jx47V6pOcnIyRI0fCyckJCoUCLVq0wCeffPLYw7MLFy7U+TN6+DCqt7c3zp8/j4MHD2p+rt7e3pptpusw4+HDhxESEgJra2tYWFigc+fO2Llzp87l7N+/H2+++SYcHR3h4OCAQYMG4datW5XaPlQ3MMyoXnjhhRcghEBMTMxj+6ampmLEiBEYOXIkduzYgd69e2P27Nn47rvvAABWVlYYO3YsNmzYUO7QZWRkJExNTTVhER8fj5CQEGRmZmLdunX48ssvcfr0abz//vs6l33q1CnMnDkTkydPxu7du/HSSy/hwoUL6Ny5M86fP4/ly5dj69at8PPzw5gxY7BkyZJHrsuRI0cwbNgw+Pr64ocffsDOnTsxf/58lJaWavrcuXMHnTt3xp49e/Dee+9hx44d6NmzJ2bMmIGJEyc+dntVxrZt2+Dr64u2bdviyJEjOHLkCLZt21Zh/4MHD6JHjx7IysrCmjVr8P3338Pa2hr9+/fXuXc9fvx4mJiYYOPGjViyZAkOHDiAkSNHVkvtJBGCSALu3LkjAIgFCxbo/Dw4OFi0bNmywul/++03AUAsXrz4kcsJDg4WAMSxY8e02v38/ERYWJjm/ZUrV4SRkZH49NNPNW0FBQXCwcFBvPrqq5q2YcOGCXNzc5GamqppKy0tFc2bNxcARGJioqbdy8tLyOVyceHCBa1lv/zyy0KhUIjk5GSt9t69ewsLCwuRmZlZ4fp8/PHHAsAj+8yaNUvnOr/55ptCJpNp1fPwz2DBggVC19dIVFRUufVr2bKlCA4OLtc3MTFRABBRUVGatk6dOgknJyeRk5OjaSstLRX+/v7Cw8NDqNVqreW89dZbWvNcsmSJACBSUlIqXG+qW7hnRvWCqMJj+1xcXNChQwettoCAAK1DlL6+vujXrx8iIyM18964cSMyMjK09mb279+PkJAQODs7a9rkcjmGDRumc9kBAQFo2rSpVtu+ffsQEhICT09PrfYxY8YgPz8fR44cqXBdnnnmGQDA0KFD8eOPP+LmzZvl+uzbtw9+fn7l1nnMmDEQQmDfvn0Vzl8f8vLycOzYMQwePBhWVlaadrlcjvDwcNy4cQMXLlzQmuaFF17Qeh8QEAAAlT6sTNLHMKN6oexLzc3N7bF9HRwcyrUpFAoUFBRotU2ZMgWXLl1CdHQ0AGDFihUICgpCu3btNH0yMjLg4uJSbn662oD7I/selpGRobO9bF0yMjIqXJfnnnsO27dvR2lpKUaNGgUPDw/4+/vj+++/r5b568O9e/cghKhSTQ//zBQKBQCU+5lR3cUwo3phx44dkMlkeO6556ptnj169IC/vz+++OILxMbG4tSpU/jXv/6l1cfBwQGpqanlptXVBkDnYAoHBwekpKSUay8b4ODo6PjIOgcMGIC9e/ciKysLBw4cgIeHB4YPH67Zo3ua+ZuZmQEAioqKtNrT09MfWdOj2NnZwcjI6KnWmeofhhnVeVFRUfjtt9/wyiuvoGHDhtU678mTJ2Pnzp2YPXs2nJ2dMWTIEK3Pu3fvjr179+L27duaNpVKVaVLBEJCQrBv375yo/PWr18PCwuLSg9rVygUCA4OxuLFiwEAp0+f1sw/Pj4ep06dKjd/mUyG7t27VzjPshGJZ86c0Wr/5ZdfdC6/MntKlpaW6NixI7Zu3arVX61W47vvvoOHh0e5Q7FEvM6MarXffvsNeXl5yMnJAXB/dOBPP/0EAOjTpw8sLCw0fQsKCnD06FHN/1+9ehXbt2/Hr7/+iuDgYHz55ZfVXt/IkSMxe/ZsHDp0CPPmzYOpqanW5/PmzcOOHTvQo0cPzJ8/HxYWFlixYgXy8vIqvYwFCxbg119/Rffu3TF//nzY29tjw4YN2LlzJ5YsWQKlUlnhtPPnz8eNGzcQEhICDw8PZGZm4rPPPoOJiQmCg4MBANOmTcP69evRt29fLFq0CF5eXti5cyciIyPx5ptvPjI4+vTpA3t7e4wbNw6LFi2CsbEx1q1bh+vXr5fr26pVK/zwww/YtGkTfH19YWZmhlatWumcb0REBHr16oXu3btjxowZMDU1RWRkJM6dO4fvv/++Wi7ZoDrGsONPiB7Ny8tLAND5enCkXNkoxLKXpaWl8PX1FYMHDxabN28WKpWqUsuraFTk6NGjhZeXl85pxowZI4yNjcWNGzd0fv7nn3+KTp06CYVCIVxcXMTMmTPFqlWrdI5m7Nu3r855nD17VvTv318olUphamoqWrdurTX6ryK//vqr6N27t3B3dxempqbCyclJ9OnTR8TExGj1u3btmhg+fLhwcHAQJiYmolmzZuKjjz4qt92gY0Tp8ePHRefOnYWlpaVwd3cXCxYsEF9//XW59UtKShKhoaHC2tpaANBsT12jGYUQIiYmRvTo0UNYWloKc3Nz0alTJ/HLL79o9SkbzfjXX39pte/fv18AEPv373/sNqK6QSZEFYZ5EZGW4uJieHt749lnn8WPP/5o6HKI6i0eZiR6Anfu3MGFCxcQFRWF27dvY9asWYYuiaheY5gRPYGdO3fi1VdfhaurKyIjI7WG4xNRzeNhRiIikjy9Ds0/dOgQ+vfvDzc3N8hkMmzfvv2x0xw8eBCBgYEwMzODr6+vXkagERFR3aLXMMvLy0Pr1q3xxRdfVKp/YmIi+vTpg65du+L06dOYM2cOJk+ejC1btuizTCIikrgaO8wok8mwbdu2Rz5k8Z133sGOHTuQkJCgaZswYQL+/vvvR95/joiI6rdaNQDkyJEjCA0N1WoLCwvDmjVrUFJSAhMTk3LTFBUVad1KR61W4+7du3BwcOCFlUREEiaEQE5ODtzc3GBk9OgDibUqzFJTU7XuLg4Azs7OKC0tRXp6us4bj0ZERGgeYEhERHXP9evX4eHh8cg+tSrMgPI3Wi07ClrRXtbs2bMxffp0zfusrCw0bNgQ169fh42Njf4KJSIivcrOzoanpyesra0f27dWhZmLi0u5u4mnpaXB2NhY52M5gPs3Ly173MODbGxsGGZERHVAZU4Z1aq75gcFBWmeDVVmz549aN++vc7zZURERICewyw3NxdxcXGIi4sDcH/ofVxcHJKTkwHcP0Q4atQoTf8JEybg2rVrmD59OhISErB27VqsWbMGM2bM0GeZREQkcXo9zHjixAmtZyGVndsaPXo01q1bh5SUFE2wAYCPjw927dqFadOmYcWKFXBzc8Py5cvx0ksv6bNMIiKSuDp3O6vs7GwolUpkZWXxnBnpnUqlQklJiaHLIJIsuVwOY2NjnefFqvJ9XqsGgBBJSW5uLm7cuIE69vcgUY2zsLCAq6truYfbVgXDjOgJqFQq3LhxAxYWFmjQoAEv0Cd6AkIIFBcX486dO0hMTESTJk0ee3F0RRhmRE+gpKQEQgg0aNAA5ubmhi6HSLLMzc1hYmKCa9euobi4GGZmZk80n1o1NJ9IarhHRvT0nnRvTGse1VAHERGRQTHMiIhI8hhmREQkeQwzIqJq9OKLL8LOzg6DBw/WtF2/fh3dunWDn58fAgICsHnzZgNWWL1ycnLwzDPPoE2bNmjVqhVWr15tkDo4mpGIqBpNnjwZY8eOxTfffKNpMzY2xrJly9CmTRukpaWhXbt26NOnDywtLQ1YafWwsLDAwYMHYWFhgfz8fPj7+2PQoEEV3hxeX7hnRkQG061bN0ydOlUv887IyICTkxOSkpL0Mv+KdO/evdwjS1xdXdGmTRsAgJOTE+zt7XH37t0K56HP7VLd5HI5LCwsAACFhYVQqVSaGwkMHjwYS5curZE6GGZE9cyYMWMgk8kgk8lgYmICZ2dn9OrVC2vXroVara7RWrZu3Yr33ntP8746v8QjIiLQv39/eHt7AwAOHTqE/v37w83NDTKZDNu3b9c5XWRkJHx8fGBmZobAwEDExMRUSz1lTpw4AbVaDU9Pz2qdryFlZmaidevW8PDwwNtvvw1HR0cAwPz58/HBBx8gOztb7zUwzIjqoeeffx4pKSlISkrCb7/9hu7du2PKlCno168fSktLa6wOe3v7Sj14saoKCgqwZs0ajB8/XtOWl5eH1q1b44svvqhwuk2bNmHq1KmYO3cuTp8+ja5du6J3795aN0QPDAyEv79/udetW7ceW1dGRgZGjRqFVatWPd0K1qDKrK+trS3+/vtvJCYmYuPGjbh9+zYAICAgAN7e3tiwYYP+CxV1TFZWlgAgsrKyDF0K1WEFBQUiPj5eFBQUGLqUKhs9erQYMGBAufa9e/cKAGL16tVCCCHUarVYvHix8PHxEWZmZiIgIEBs3rxZa5rg4GAxadIkMXPmTGFnZyecnZ3FggULtPps3rxZ+Pv7CzMzM2Fvby9CQkJEbm6uZvopU6Zo6gKg9Xr33XeFvb29KCws1JrnoEGDRHh4eIXruGXLFuHo6Fjh5wDEtm3byrV36NBBTJgwQautefPmYtasWRXOS5f9+/eLl156SautsLBQdO3aVaxfv/6x0z+4XYQQ4rfffhM2Njbim2++EUIIkZ2dLYYPHy4sLCyEi4uLWLp0ablpKprvxIkTxZQpU4Stra1wcnISX331lcjNzRVjxowRVlZWwtfXV+zatatK6/ugCRMmiB9//FHzfuHChaJr166PnKaif09V+T7nnhkRAQB69OiB1q1bY+vWrQCAefPmISoqCitXrsT58+cxbdo0jBw5EgcPHtSa7ptvvoGlpSWOHTuGJUuWYNGiRZqH7KakpOCVV17B2LFjkZCQgAMHDmDQoEE6b8782WefISgoCK+99hpSUlKQkpKCf//731CpVNixY4emX3p6On799Ve8+uqrFa7LoUOH0L59+yqtf3FxMU6ePInQ0FCt9tDQUMTGxlZpXg8TQmDMmDHo0aMHwsPDqzTtDz/8gKFDh2L9+vWa5z9Onz4df/75J3bs2IHo6GjExMTg1KlTlZrfN998A0dHRxw/fhyTJk3Cm2++iSFDhqBz5844deoUwsLCEB4ejvz8/ErN7/bt25rDiNnZ2Th06BCaNWum+bxDhw44fvw4ioqKqrTeVcXRjETVQAiBghKVQZZtbiKvtttqNW/eHGfOnEFeXh6WLl2Kffv2ISgoCADg6+uLw4cP46uvvkJwcLBmmoCAACxYsAAA0KRJE3zxxRfYu3cvevXqhZSUFJSWlmLQoEHw8vICALRq1UrnspVKJUxNTWFhYQEXFxdN+/DhwxEVFYUhQ4YAADZs2AAPDw9069atwvVISkqCm5tbldY9PT0dKpUKzs7OWu3Ozs5ITU2t9HzCwsJw6tQp5OXlwcPDA9u2bUNRURE2bdqEgIAAzbm6b7/9tsJtUSYyMhJz5szBzz//rHk2ZE5ODr755hts3LgRISEhAICoqKhKr2/r1q0xb948APcfkPzhhx/C0dERr732GoD757lWrlyJM2fOoFOnTo+d340bNzBu3DgIISCEwMSJExEQEKD53N3dHUVFRUhNTdX8DugDw4yoGhSUqOA3/3eDLDt+URgsTKvnn7IQAjKZDPHx8SgsLESvXr20Pi8uLkbbtm212h784gLuj9xLS0sDcP+LMyQkBK1atUJYWBhCQ0MxePBg2NnZVbqm1157Dc888wxu3rwJd3d3REVFaQaxVKSgoOCJb1j78HzLtkll/f677t+Dqg6u2bJlC27fvo3Dhw+jQ4cOmvarV6+ipKREq02pVGrtDT3Kgz8vuVwOBwcHrVAtC/Oyn+HjBAYGIi4ursLPy27EXdk9vSfFw4xEpJGQkAAfHx/NF+/OnTsRFxenecXHx+Onn37SmsbExETrvUwm00wvl8sRHR2N3377DX5+fvj888/RrFkzJCYmVrqmtm3bonXr1li/fj1OnTqFs2fPYsyYMY+cxtHREffu3av0Msqmkcvl5fbC0tLSyu2t1YQ2bdqgQYMGiIqK0josW/b/ukK3MnT9vB5sK5tvdY1sLbsEoUGDBtUyv4pwz4yoGpibyBG/KMxgy64O+/btw9mzZzFt2jT4+flBoVAgOTlZ65Dik5DJZOjSpQu6dOmC+fPnw8vLC9u2bcP06dPL9TU1NYVKVf5w7fjx4/Hpp5/i5s2b6Nmz52OHtbdt2xbfffddleo0NTVFYGAgoqOj8eKLL2rao6OjMWDAgCrNqzo0atQIn3zyCbp16wa5XK4ZhdmoUSOYmJjg+PHjmu2QnZ2NS5cuPfXPSh/OnTsHDw8PzXB9fWGYEVUDmUxWbYf6akLZOQyVSoXbt29j9+7diIiIQL9+/TBq1CjI5XLMmDED06ZNg1qtxrPPPovs7GzExsbCysoKo0ePrtRyjh07hr179yI0NBROTk44duwY7ty5gxYtWujs7+3tjWPHjiEpKQlWVlawt7eHkZERRowYgRkzZmD16tVYv379Y5cbFhaG2bNn4969e5pDmrm5ubh8+bKmT2JiIuLi4mBvb4+GDRsCuD+wIjw8HO3bt0dQUBBWrVqF5ORkTJgwoVLrW92aNm2K/fv3o1u3bpq7iFhbW2P06NGYOXMm7O3t4eTkhAULFsDIyKhWPpIoJiam3KAafZDOvz4iqja7d++Gq6srjI2NYWdnh9atW2P58uUYPXq05tlS7733HpycnBAREYGrV6/C1tYW7dq1w5w5cyq9HBsbGxw6dAjLli1DdnY2vLy88Mknn6B37946+8+YMQOjR4+Gn58fCgoKkJiYCG9vb9jY2OCll17Czp07MXDgwMcut1WrVmjfvj1+/PFHvPHGGwDuX6xcNogCgGbPcPTo0Vi3bh0AYNiwYcjIyMCiRYuQkpICf39/7Nq1S68DFx6nWbNm2Ldvn2YP7ZNPPsHSpUsxYcIE9OvXDzY2Nnj77bdx/fr1Jz5PqC+FhYXYtm1bhecRq5NMVPZAq0RkZ2dDqVQiKysLNjY2hi6H6qjCwkIkJiZq7hRB+terVy+0aNECy5cvr1T/Xbt2YcaMGTh37ly1PPyxNsvLy4O7uzs++eQTjBs3ztDlaKxYsQI///wz9uzZ88h+Ff17qsr3OffMiKhWu3v3Lvbs2YN9+/Y98u4dD+vTpw8uXbqEmzdv1qlbRwHA6dOn8c8//6BDhw7IysrCokWLAMAg5/YexcTEBJ9//nmNLIthRkS1Wrt27XDv3j0sXry40sPPy0yZMkVPVRnexx9/jAsXLmgGrsTExOh9kEVVvf766zW2LIYZEdVqNX3Xeylo27YtTp48aegyapW6fSCZiIjqBYYZERFJHsOMiIgkj2FGRESSxzAjIiLJY5gREZHkMcyIiEjyGGZERCR5DDMiIpI8hhkREUkew4yIiCSPYUZERJLHMCOiJ9KtWzdMnTq11s6vTEZGBpycnGrkhsX6WoeqGjx4MJYuXWroMmoUw4yoHunfvz969uyp87MjR45AJpPh1KlTNVzVfVu3bsV7770HoHpDISIiAv3794e3t7em7dChQ+jfvz/c3Nwgk8mwfft2ndNGRkZqHhhZ9pgVQxozZgxmzZr12H7z58/HBx98gOzs7BqoqnbQe5hV9Zdhw4YNaN26NSwsLODq6opXX30VGRkZ+i6TqF4YN24c9u3bh2vXrpX7bO3atWjTpg3atWtngMoAe3t7WFtbV+s8CwoKsGbNGowfP16rPS8vD61bt37kwz43bdqEqVOnYu7cuTh9+jS6du2K3r17Izk5uVprrCy1Wo2dO3dW6gGcAQEB8Pb2xoYNG2qgslpC6NEPP/wgTExMxOrVq0V8fLyYMmWKsLS0FNeuXdPZPyYmRhgZGYnPPvtMXL16VcTExIiWLVuKgQMHVnqZWVlZAoDIysqqrtUgKqegoEDEx8eLgoICQ5dSJSUlJcLZ2VksXLhQqz0vL09YW1uLzz//XAghhFqtFosXLxY+Pj7CzMxMBAQEiM2bN2tNExwcLKZMmSKEEKKwsFBMmjRJNGjQQCgUCtGlSxdx/Phxrf4qlUp8+OGHolGjRsLU1FR4enqK999/v9z8Ro8eLQBovRITE8U333wj7O3tRWFhodZ8Bw0aJMLDw3Wu75YtW4Sjo+MjtwkAsW3btnLtHTp0EBMmTNBqa968uZg1a1aF83pwmwghxG+//SZsbGzEN998o/l84sSJYsqUKcLW1lY4OTmJr776SuTm5ooxY8YIKysr4evrK3bt2lVu3ocOHRJOTk5CpVKJzZs3C39/f2FmZibs7e1FSEiIyM3N1eq/cOFC0bVr10eue21R0b+nqnyf6zXMqvrL8NFHHwlfX1+ttuXLlwsPD49KL5NhRjWh3D8+tVqIolzDvNTqKtU+c+ZM4e3tLdQPTLdu3TqhUCjE3bt3hRBCzJkzRzRv3lzs3r1bXLlyRURFRQmFQiEOHDigmebBL+7JkycLNzc3sWvXLnH+/HkxevRoYWdnJzIyMjT93377bWFnZyfWrVsnLl++LGJiYsTq1avLzS8zM1MEBQWJ1157TaSkpIiUlBRRWloq8vPzhVKpFD/++KNmmjt37ghTU1Oxb98+nes6ZcoU8fzzzz9ye+gKs6KiIiGXy8XWrVu12idPniyee+65Cuf14Db5/vvvhbW1tdi+fbvW59bW1uK9994TFy9eFO+9954wMjISvXv3FqtWrRIXL14Ub775pnBwcBB5eXla854xY4YYN26cuHXrljA2NhZLly4ViYmJ4syZM2LFihUiJydHq/+uXbuEQqEoF/61UXWEmd6eNF1cXIyTJ0+WO74bGhqK2NhYndN07twZc+fOxa5du9C7d2+kpaXhp59+Qt++fStcTlFREYqKijTv69MxYqpFSvKB/7gZZtlzbgGmlpXuPnbsWHz00Uc4cOAAunfvDuD+IcZBgwbBzs4OeXl5WLp0Kfbt24egoCAAgK+vLw4fPoyvvvoKwcHBWvPLy8vDypUrsW7dOvTu3RsAsHr1akRHR2PNmjWYOXMmcnJy8Nlnn+GLL77A6NGjAQCNGjXCs88+W64+pVIJU1NTWFhYwMXFRdNubm6O4cOHIyoqCkOGDAFw/7SEh4cHunXrpnNdk5KS4OZW9Z9Leno6VCoVnJ2dtdqdnZ2Rmpr62OkjIyMxZ84c/Pzzz5ptXKZ169aYN28eAGD27Nn48MMP4ejoiNdeew3A/fNdK1euxJkzZ9CpUyfNdDt27MDHH3+MlJQUlJaWYtCgQfDy8gIAtGrVqlwN7u7uKCoqQmpqqqZfXaa3MHuSX4bOnTtjw4YNGDZsGAoLC1FaWooXXngBn3/+eYXLiYiIwLvvvluttRPVZc2bN0fnzp2xdu1adO/eHVeuXEFMTAz27NkDAIiPj0dhYSF69eqlNV1xcTHatm1bbn5XrlxBSUkJunTpomkzMTFBhw4dkJCQAABISEhAUVERQkJCnqr21157Dc888wxu3rwJd3d3REVFYcyYMZDJZDr7FxQUwMzM7ImX9/B8hRAVLqvMli1bcPv2bRw+fBgdOnQo93lAQIDm/+VyORwcHLTCqOw7My0tTdOWkJCAGzduoGfPnjA1NUVISAhatWqFsLAwhIaGYvDgwbCzs9Najrm5OQAgPz+/kmsrbXoLszJV+WWIj4/H5MmTMX/+fISFhSElJQUzZ87EhAkTsGbNGp3TzJ49G9OnT9e8z87OhqenZ/WtAFFlmFjc30My1LKraNy4cZg4cSJWrFiBqKgoeHl5aYJGrVYDAHbu3Al3d3et6RQKRbl5CSEAPPrfetkX69Nq27YtWrdujfXr1yMsLAxnz57FL7/8UmF/R0dH3Lt3r8rLcXR0hFwuL/eHd1paWrk/0B/Wpk0bnDp1ClFRUXjmmWfKbRcTExOt9zKZTKutrH/ZzwG4v1fWq1cvzXaMjo5GbGws9uzZg88//xxz587FsWPH4OPjo5nm7t27AIAGDRpUdrUlTW+jGZ/klyEiIgJdunTBzJkzERAQgLCwMERGRmLt2rVISUnROY1CoYCNjY3Wi6jGyWT3D/UZ4vWYPQVdhg4dCrlcjo0bN+Kbb77Bq6++qvkS9fPzg0KhQHJyMho3bqz10vWHYuPGjWFqaorDhw9r2kpKSnDixAm0aNECANCkSROYm5tj7969larP1NQUKpVK52fjx49HVFQU1q5di549ez7yj9e2bdsiPj6+Ust8ePmBgYGIjo7Wao+Ojkbnzp0fOW2jRo2wf/9+/Pzzz5g0aVKVl63Lzz//jBdeeEHzXiaToUuXLnj33Xdx+vRpmJqaYtu2bVrTnDt3Dh4eHnB0dKyWGmo7ve2ZPfjL8OKLL2rao6OjKxxamp+fD2Nj7ZLkcjmA//31R0RPz8rKCsOGDcOcOXOQlZWFMWPGaD6ztrbGjBkzMG3aNKjVajz77LPIzs5GbGwsrKysNOe8ylhaWuLNN9/EzJkzYW9vj4YNG2LJkiXIz8/HuHHjAABmZmZ455138Pbbb8PU1BRdunTBnTt3cP78eU2fB3l7e+PYsWNISkqClZUV7O3tYWR0/2/vESNGYMaMGVi9ejXWr1//yPUMCwvD7Nmzce/ePa3DcLm5ubh8+bLmfWJiIuLi4jT1A8D06dMRHh6O9u3bIygoCKtWrUJycjImTJjw2O3btGlT7N+/H926dYOxsTGWLVv22GkqkpaWhr/++ktzLdyxY8ewd+9ehIaGwsnJCceOHcOdO3c0fziUiYmJQWho6BMvV3L0MDBFo2xo/po1a0R8fLyYOnWqsLS0FElJSUIIIWbNmqU1pDYqKkoYGxuLyMhIceXKFXH48GHRvn170aFDh0ovk6MZqSZIdWj+g2JjYwUAERoaWu4ztVotPvvsM9GsWTNhYmIiGjRoIMLCwsTBgwc1fR4cuVdQUCAmTZokHB0dHzk0//333xdeXl7CxMRENGzYUPznP//ROb8LFy6ITp06CXNzc83Q/AeFh4frHKavS6dOncSXX36p1bZ///5yw/8BiNGjR2v1W7FihfDy8hKmpqaiXbt2Wuuvy8ND8+Pj44WTk5OYPn26zs+FEMLLy0t8+umnWm14YITl119/Lbp06aI1z7CwMM1lEE2bNtVcUlGmoKBA2NjYiCNHjjyy3tqi1g/NF+LRvwyjR48WwcHBWv2XL18u/Pz8hLm5uXB1dRUjRowQN27cqPTyGGZUE+pCmElZz549xaRJkyrVd+fOnaJFixZCpVLpuSr96N+/v1i8eHGVpvniiy9Er1699FRR9auOMJMJUbeO32VnZ0OpVCIrK4vnz0hvCgsLkZiYqLm7DdWMu3fvYs+ePRgxYgTi4+PRrFmzSk332WefYdCgQZIcHLZkyRK88sorVap91apVCA4OrvT2MbSK/j1V5ftc76MZiYiqS7t27XDv3j0sXry4Sl/UU6ZM0WNV+vX2229XeZrXX39dD5XUbgwzIpKMmrjzPUkT75pPRESSxzAjIiLJY5gREZHkMcyIiEjyGGZERCR5DDMiIpI8hhkREUkew4yIiCSPYUZERJLHMCOiemXNmjV1+tEogwcPxtKlSw1dRo1jmBHVI2PGjIFMJtM83djX1xczZsxAXl6eoUurEUVFRZg/fz7+7//+T6t9y5YtmoeS+vn5lXvQpS5nz55FcHAwzM3N4e7ujkWLFmk9d/Hw4cPo0qULHBwcYG5ujubNm+PTTz995DyTkpI0Px+ZTAY7Ozs899xzOHjwYKXXcf78+fjggw+QnZ1d6WnqAoYZUT3z/PPPIyUlBVevXsX777+PyMhIzJgxQ2ffkpKSGq7u8Z6mpi1btsDKygpdu3bVtB05cgTDhg1DeHg4/v77b4SHh2Po0KE4duxYhfPJzs5Gr1694Obmhr/++guff/45Pv74Y609IktLS0ycOBGHDh1CQkIC5s2bh3nz5mHVqlWPrfOPP/5ASkoKDh48CBsbG/Tp0weJiYmVWseAgAB4e3tjw4YNlepfZ+jl4TQGxOeZUU14+PlLarVa5BXnGeSlVqsrXffo0aPFgAEDtNrGjx8vXFxchBBCLFiwQLRu3VqsWbNG+Pj4CJlMJtRqtcjMzBSvvfaaaNCggbC2thbdu3cXcXFxmnnExcWJbt26CSsrK2FtbS3atWsn/vrrLyGEEElJSaJfv37C1tZWWFhYCD8/P7Fz504hxP0H8iqVSq16tm3bJh78anrSmnTp37+/mDFjhlbb0KFDxfPPP6/VFhYWJl5++eUK5xMZGSmUSqXWw0EjIiKEm5vbI38eL774ohg5cmSFnycmJgoA4vTp05q2GzduCADiyy+/FOnp6eLll18W7u7uwtzcXPj7+4uNGzeWm8/ChQtF165dK1xObVMdzzPjXfOJqkFBaQE6buxokGUfG34MFiYWTzy9ubm51t7O5cuX8eOPP2LLli2Qy+UAgL59+8Le3h67du2CUqnEV199hZCQEFy8eBH29vYYMWIE2rZti5UrV0IulyMuLg4mJiYAgH/9618oLi7GoUOHYGlpifj4eFhZWVWpxiepSZeYmBiMGDFCq+3IkSOYNm2aVltYWBiWLVtWYT1HjhxBcHAwFAqF1jSzZ89GUlISfHx8yk1z+vRpxMbG4v3336/sagMALCzu/2xLSkpQWFiIwMBAvPPOO7CxscHOnTsRHh4OX19fdOz4v9+/Dh06ICIiAkVFRVo11mUMM6J67Pjx49i4cSNCQkI0bcXFxfj222/RoEEDAMC+fftw9uxZpKWlab4YP/74Y2zfvh0//fQTXn/9dSQnJ2PmzJlo3rw5AKBJkyaa+SUnJ+Oll15Cq1atAAC+vr5VrvNJanpYZmYmMjMz4ebmptWempoKZ2dnrTZnZ2ekpqZWWE9qaiq8vb3LTVP22YNh5uHhgTt37qC0tBQLFy7E+PHjK73eeXl5mD17NuRyOYKDg+Hu7q51SHjSpEnYvXs3Nm/erBVm7u7uKCoqQmpqKry8vCq9PCljmBFVA3NjcxwbXvE5Fn0vuyp+/fVXWFlZobS0FCUlJRgwYAA+//xzzedeXl6a0ACAkydPIjc3Fw4ODlrzKSgowJUrVwAA06dPx/jx4/Htt9+iZ8+eGDJkCBo1agQAmDx5Mt58803s2bMHPXv2xEsvvYSAgIAq1fwkNT2soKAAAHQ+GVwmk2m9F0KUa6vMNLraY2JikJubi6NHj2LWrFlo3LgxXnnllUfOu3PnzjAyMkJ+fj5cXV2xbt06tGrVCiqVCh9++CE2bdqEmzdvoqioCEVFRbC0tNSa3tz8/u9Efn7+I5dTlzDMiKqBTCZ7qkN9Nal79+5YuXIlTExM4ObmpjkcWObhL0a1Wg1XV1ccOHCg3LxsbW0BAAsXLsTw4cOxc+dO/Pbbb1iwYAF++OEHvPjiixg/fjzCwsKwc+dO7NmzBxEREfjkk08wadIkGBkZaY0ABHQP8HiSmh7m4OAAmUyGe/fuabW7uLiU2wtLS0srt7dWmWkAlJuubC+tVatWuH37NhYuXPjYMNu0aRP8/Pxga2urFdiffPIJPv30UyxbtgytWrWCpaUlpk6diuLiYq3p7969CwBafwDUdRzNSFTPWFpaonHjxvDy8ioXZLq0a9cOqampMDY2RuPGjbVejo6Omn5NmzbFtGnTsGfPHgwaNAhRUVGazzw9PTFhwgRs3boV//73v7F69WoA979sc3JytC4NiIuLq7aaHmRqago/Pz/Ex8drtQcFBSE6Olqrbc+ePejcuXOFyw8KCsKhQ4e0QmTPnj1wc3Mrd/jxQUIIFBUVPXb9PD090ahRo3J7njExMRgwYABGjhyJ1q1bw9fXF5cuXSo3/blz5+Dh4VHhtqiLGGZE9Eg9e/ZEUFAQBg4ciN9//x1JSUmIjY3FvHnzcOLECRQUFGDixIk4cOAArl27hj///BN//fUXWrRoAQCYOnUqfv/9dyQmJuLUqVPYt2+f5rOOHTvCwsICc+bMweXLl7Fx40asW7fuqWuqSFhYGA4fPqzVNmXKFOzZsweLFy/GP//8g8WLF+OPP/7A1KlTNX2++OILrfOKw4cPh0KhwJgxY3Du3Dls27YN//nPfzB9+nTNYcYVK1bgl19+waVLl3Dp0iVERUXh448/xsiRIyu76ctp3LgxoqOjERsbi4SEBLzxxhs6z+3FxMTU6QvDddLHMEtD4tB8qgkVDSWu7XQNzX9Q2TD4h2VnZ4tJkyYJNzc3YWJiIjw9PcWIESNEcnKyKCoqEi+//LLw9PQUpqamws3NTUycOFGzbSZOnCgaNWokFAqFaNCggQgPDxfp6emaeW/btk00btxYmJmZiX79+olVq1bpHJpflZoqkpCQIMzNzUVmZqZW++bNm0WzZs2EiYmJaN68udiyZUu57eLl5aXVdubMGdG1a1ehUCiEi4uLWLhwodaw/OXLl4uWLVsKCwsLYWNjI9q2bSsiIyOFSqWqsD5dQ/MflJGRIQYMGCCsrKyEk5OTmDdvnhg1apTWz7SgoEDY2NiII0eOVLic2qY6hubLhHjogLXEZWdnQ6lUIisrCzY2NoYuh+qowsJCJCYmwsfHR+eAAqq9hg4dirZt22L27NmGLkUvVqxYgZ9//hl79uwxdCmVVtG/p6p8n/MwIxHVKx999FGVr3OTEhMTE63RqfUFRzMSUb3i5eWFSZMmGboMvdF1jV19wD0zIiKSPIYZ0VOoY6eciQyiOv4dMcyInkDZ/QEfvliViKqu7E4llbnusSI8Z0b0BIyNjWFhYYE7d+7AxMQERkb8u5CoqoQQyM/PR1paGmxtbTV/JD4JhhnRE5DJZHB1dUViYiKuXbtm6HKIJM3W1hYuLi5PNQ+GGdETMjU1RZMmTXiokegpmJiYPNUeWRmGGdFTMDIy4kXTRLUAD/QTEZHkMcyIiEjyGGZERCR5DDMiIpI8hhkREUkew4yIiCRP72EWGRmpeUZNYGAgYmJiHtm/qKgIc+fOhZeXFxQKBRo1aoS1a9fqu0wiIpIwvV5ntmnTJkydOhWRkZHo0qULvvrqK/Tu3Rvx8fFo2LChzmmGDh2K27dvY82aNWjcuDHS0tJQWlqqzzKJiEji9Pqk6Y4dO6Jdu3ZYuXKlpq1FixYYOHAgIiIiyvXfvXs3Xn75ZVy9ehX29vZPtEw+aZqIqG6oFU+aLi4uxsmTJxEaGqrVHhoaitjYWJ3T7NixA+3bt8eSJUvg7u6Opk2bYsaMGSgoKKhwOUVFRcjOztZ6ERFR/aK3w4zp6elQqVRwdnbWand2dkZqaqrOaa5evYrDhw/DzMwM27ZtQ3p6Ot566y3cvXu3wvNmERERePfdd6u9fiIikg69DwCRyWRa74UQ5drKqNVqyGQybNiwAR06dECfPn2wdOlSrFu3rsK9s9mzZyMrK0vzun79erWvAxER1W562zNzdHSEXC4vtxeWlpZWbm+tjKurK9zd3aFUKjVtLVq0gBACN27cQJMmTcpNo1AooFAoqrd4IiKSFL3tmZmamiIwMBDR0dFa7dHR0ejcubPOabp06YJbt24hNzdX03bx4kUYGRnBw8NDX6USEZHE6fUw4/Tp0/H1119j7dq1SEhIwLRp05CcnIwJEyYAuH+IcNSoUZr+w4cPh4ODA1599VXEx8fj0KFDmDlzJsaOHQtzc3N9lkpERBKm1+vMhg0bhoyMDCxatAgpKSnw9/fHrl274OXlBQBISUlBcnKypr+VlRWio6MxadIktG/fHg4ODhg6dCjef/99fZZJREQSp9frzAyB15kREdUNteI6MyIioprCMCMiIsljmBERkeQxzIiISPIYZkREJHkMMyIikjyGGRERSR7DjIiIJI9hRkREkscwIyIiyWOYERGR5DHMiIhI8hhmREQkeQwzIiKSPIYZERFJHsOMiIgkj2FGRESSxzAjIiLJY5gREZHkMcyIiEjyGGZERCR5DDMiIpI8hhkREUkew4yIiCSPYUZERJLHMCMiIsljmBERkeQxzIiISPIYZkREJHkMMyIikjyGGRERSR7DjIiIJI9hRkREkscwIyIiyWOY6XAn7x6Wn1qOEnWJoUshIqJKYJg9RC3UGLh1JFafXY1nvx2Ao8mXDF0SERE9BsPsIUYyI9gV9wEA5OM6xv8xHG/98hkKSooNXBkREVWEYabDT6MmYKr/RzBTN4RMXoiYu1+j8/oXsCHugKFLIyIiHWRCCGHoIqpTdnY2lEolsrKyYGNj81TzKlWVYu7er7HrxlpAXgAAcDHqgs/C5sHPyaM6yiUiogpU5ftc73tmkZGR8PHxgZmZGQIDAxETE1Op6f78808YGxujTZs2+i3wEYzlxlgcOgE7XvwFDU26QwgZUtV/YujOgXhjxyfIKy40WG1ERPQ/eg2zTZs2YerUqZg7dy5Onz6Nrl27onfv3khOTn7kdFlZWRg1ahRCQkL0WV6l+dg5Y+fw5Xivw1dQqLwhMypC7L116PJtP6w8tsvQ5RER1Xt6PczYsWNHtGvXDitXrtS0tWjRAgMHDkRERESF07388sto0qQJ5HI5tm/fjri4uEovszoPM+pSqlJhwf512JH8NSDPBQDYifb4qMdcdGzYuNqXR0RUX9WKw4zFxcU4efIkQkNDtdpDQ0MRGxtb4XRRUVG4cuUKFixYUKnlFBUVITs7W+ulT8ZyOT7oOQ67B/+KpubPQwgZ7slOYNzeIXhl80Kk5ep3+UREVJ7ewiw9PR0qlQrOzs5a7c7OzkhNTdU5zaVLlzBr1ixs2LABxsbGlVpOREQElEql5uXp6fnUtVeGu40Dtgz9CJ91XQ8r0Qwyo1Kcy9+CkB/7YMG+dVCpVTVSBxER1cAAEJlMpvVeCFGuDQBUKhWGDx+Od999F02bNq30/GfPno2srCzN6/r1609dc1WENGqDP0f9iLGNF8JI5QDIs7D1+ifouG4ANp87VKO1EBHVV3o7Z1ZcXAwLCwts3rwZL774oqZ9ypQpiIuLw8GDB7X6Z2Zmws7ODnK5XNOmVqshhIBcLseePXvQo0ePxy5X3+fMHiWnqADTdn+Bo3c3QWZUBABwMuqEj3rMRjt33xqthYhI6mrFOTNTU1MEBgYiOjpaqz06OhqdO3cu19/GxgZnz55FXFyc5jVhwgQ0a9YMcXFx6Nixo75KrTbWCnN8PWAmNvXZBnd5NwghQ5r6KEbteQmvbF6I1JwsQ5dIRFQnVe7E1BOaPn06wsPD0b59ewQFBWHVqlVITk7GhAkTANw/RHjz5k2sX78eRkZG8Pf315reyckJZmZm5dpru5bOntg98nP8+s8JvHckAvlGF3Eufwt6bf4dIS4j8Z9e42FhojB0mUREdYZez5kNGzYMy5Ytw6JFi9CmTRscOnQIu3btgpeXFwAgJSXlsdecSVm/5u1xZPRmvNp4PuSqBoA8F3vvfImg9b2xJOZHqNVqQ5dIRFQn8HZWNaSgpAhz//ga0Snfaa5PM1P7YFrgNAwP6G7g6oiIap+qfJ8zzGpYak4mZvy+HHE5P0NmdP9O/LZojfldZqBX4zaGLY6IqBZhmNXiMCtz/vZ1vP3HUlwr2QeZTA0hZHCTP4sPuk3HM568kwgREcNMAmFW5lBiPBbEfIR0cQIAIIQcPqY98GHINLR0rpkLwImIaiOGmYTCrMz2+FgsOf4ZcmTxAAChNkYLy+expOcU+Ni7GLg6IqKaxzCTYJiV+ebUH1gRtwIF8sv3G9SmCLDphw9DJsHT1tGwxRER1SCGmYTDDLh/55PIY7uwNn4lSoz/e+mC2gxtbF7A4p4T4aa0M2yBREQ1gGEm8TAro1Kp8Wnsdmy4uAqlxjf/22iBtsr++CDkTXjaOhi2QCIiPWKY1ZEwK1OqUuHjPzdj0+WvUSq/DQAQajO0tumL//R4C152TgaukIio+jHM6liYlSkpLcXHsT/ix8vrUCpPAQAItSn8rXvj/e7/QmMHVwNXSERUfRhmdTTMypSqVPg0diu+vxiFEuP7j7wRamM0Me+Jd4MnIsDVy8AVEhE9PYZZHQ+zMiqVGsuP7sB3/6xBsXESAECo5Who2g1zn30TXbyaGbZAIqKnwDCrJ2FWRq1WY/WJ37H2/GrkG10CAAghQwOjZzC1/esY4Ff7H59DRPQwhlk9C7MHbfx7PyLjViMLZzVtVuoWGOs/BuMCn4eRkd4fLk5EVC0YZvU4zMr8fvEUPjm2CrdURyCT3X/UjKnKAwN9R2BGlyEw5/PUiKiWY5gxzDRO3LiC9w9/hcsFezV36ZeV2qOr0yD8X7cxcLFWGrhCIiLdGGYMs3IS76bh3YNrcPLeDs3z1ITKHM0te2J2l/EI9PA1cIVERNoYZgyzCmUW5OE/Md9hz41NUMnvAACEMEIDo0C80WY0hrV6DjKZzMBVEhExzBhmlVCqKsWXf+3EhoQNyDVK0LQrVA3R32cYZnQeBksFz6sRkeEwzBhmVfL7pdNYdmwtrpcchsyo9H6jyhptbftgTtexaN7AzbAFElG9xDBjmD2RS+mp+E9MFE7c2wnIswDcvwjbWd4R41uPwLBWz3JoPxHVGIYZw+yp5JcUYemfP+HnxB9QaJSkaTcudUdPjxfxzrOvwNGS25aI9IthxjCrNjsSjuHLU98iufjP/x2CVCvQyDwYk54ZhZBGrQ1bIBHVWQwzhlm1u56ZgQ8Pf4vDt3+B2jhN026uaoK+3oMwrfNLsDEzN2CFRFTXMMwYZnqjUqmx+sQebPxnE+6KU5q7i0BlhSYWwZj4zAj0aNTKsEUSUZ3AMGOY1Yhzqcn4+Mh6nLq3G+K/A0YAwEzli16e/fHvzsPgYGltwAqJSMoYZgyzGlVUWoLVJ3Zh84UtyBB//29vTa2Al6IzXm09FC+26MyRkERUJQwzhpnBJKTdwNIjG3E8fTfUxnc07cYqF3R0fB7Tgl5BM163RkSVwDBjmBmcSqXG+tP7sTFhM1JKj0NmVALg/q2z7BCAPj798FbH/lCaWRi4UiKqrRhmDLNa5Vb2XXx6ZDP23/wVRfKk/32gMoeXWRBG+g/CUP+uPAxJRFoYZgyzWuvA1TP46tRmnM8+ACHP1LQbqRzQxjYEbwQORmevFoYrkIhqDYYZw6zWKyktxfq4vdiUsB23So9rnrUGAAqVD7q4hGFih5fQxNHFgFUSkSExzBhmkpKWm43Pj27H3uu/IVt2HjLZ/V9JIYxgI/zQwzMUb3UYADcbewNXSkQ1iWHGMJOshLQbWHF8M46kRaNYfl3TLtRyOBgFoLfP83jjmRdgZ25lwCqJqCYwzBhmdcKf1+Kx+tQ2xN09AJVxqqZdqE3gbNwWz3v3xuvP9OGISKI6imHGMKtT1Go19lz6G+vObEN8dgyEcbrmM6FWwEUeiOd9wjC+/fOwNWewEdUVDDOGWZ2lUqmxPeEYvj//Cy7mxWiNiBRqBRrIW6OXVwjGB/aFk5XScIUS0VNjmDHM6gWVWoUt8bHYFP8LLuUeeSjYjGFv1BLdPEIwrl1feNk5Ga5QInoiDDOGWb2jVqux45/j+OH8TiRk/6l1Ky0hjGAlmqKjUzBebdsPbdy8DVcoEVVaVb7P9X7LhcjISPj4+MDMzAyBgYGIiYmpsO/WrVvRq1cvNGjQADY2NggKCsLvv/+u7xKpDjAyMsJAv074Ych7OP3qH/goaD3aWA2DcakHZDI18oz+wb70rxAe3R/t1vTDqK0f4Nd/TkCtVhu6dCKqBnrdM9u0aRPCw8MRGRmJLl264KuvvsLXX3+N+Ph4NGzYsFz/qVOnws3NDd27d4etrS2ioqLw8ccf49ixY2jbtm2llsk9M3rYseSL+ObvX3Ei/RAKjK5ofSZT2cLbvD16N+qBEQEhsOHISKJao9YcZuzYsSPatWuHlStXatpatGiBgQMHIiIiolLzaNmyJYYNG4b58+dXqj/DjB7lUvpNRJ3+DYdvHcRd9XnNDZABAGoT2Bv5o7NrV4S36Q0/Jw/DFUpEVfo+N9ZXEcXFxTh58iRmzZql1R4aGorY2NhKzUOtViMnJwf29hXf+aGoqAhFRUWa99nZ2U9WMNULTRzd8Z9e4wGMR2ZBHr6N24vdV/chufAkYJyJuziNX1NO49eU5TApbQg/2054oUkPDGjRAQoTE0OXT0QV0FuYpaenQ6VSwdnZWavd2dkZqampFUyl7ZNPPkFeXh6GDh1aYZ+IiAi8++67T1Ur1U+25paYFPQCJgW9AJVKjZ0XTmLLP9E4n3kERfIklBgn4+/cZPx9+ke8d8ICTiYB6OzWGSMDQtGsgbuhyyeiB+gtzMrIZDKt90KIcm26fP/991i4cCF+/vlnODlVPKx69uzZmD59uuZ9dnY2PD09n7xgqpfkciO84PcMXvB7BgBwMf0mvv17N/68eRh3VOcAeT7S1Eex/cZRbL+xFKYqdzS1aY9Q32AMadkVVgozA68BUf2mtzBzdHSEXC4vtxeWlpZWbm/tYZs2bcK4ceOwefNm9OzZ85F9FQoFFArFU9dL9KCmju54L2QcgHEoKCnC9vgj+PXSAfyT9ReK5NdRLL+Jc3k3ce7sz/jkb1PYylqgbYOOeKlFCJ7zbs5nsxHVML0PAAkMDERkZKSmzc/PDwMGDKhwAMj333+PsWPH4vvvv8fAgQOrvEwOACF9u3AnBd/9HY2jKbFILTkDyHO0PpeV2sFV4Y9Orp0wzL87/Jx5pIDoSdSa0YxlQ/O//PJLBAUFYdWqVVi9ejXOnz8PLy8vzJ49Gzdv3sT69esB3A+yUaNG4bPPPsOgQYM08zE3N4dSWblbEzHMqCaVqlTYc/k0fr6wH2cyjiFHdhkymUqrj7HKGV4WrfGcRxCG+HeDp62jgaolkpZaE2bA/YumlyxZgpSUFPj7++PTTz/Fc889BwAYM2YMkpKScODAAQBAt27dcPDgwXLzGD16NNatW1ep5THMyJAyC3Kx+VwM9ib9iUs5cSgyStY8nw0AhJBBofaAr2UbPNcwCEP8g+Fizd9TIl1qVZjVNIYZ1SbX7qXjh7MHcPjGEVwvOKP1KBvg/q22zFRe8LFuha6eHfBSy65w50NIiQAwzBhmVGvFp93AT+cP4GjKMdwsPAO1/K7W5/f33NzhZemPILdn8KJfVzR2cDVQtUSGxTBjmJFExKVcwfaEwzieegK3CuKhMk4r10de6gQ3s5Zo59wO/Zs+iw6ejSp1eQuR1DHMGGYkUfG3b2BbQgyO3jqBGwXnUSK/pXXODQBQagsH4yZo6RCA7l7PoHfT9rDk5SlUBzHMGGZUR1zPTMe2+D9x+OZxJOaeRYHsGmQy7Tv9C7UxLOAFHys/dHRth77NgniHEqoTGGYMM6qjMgtz8es/x3Dw2gn8k3kOmaqLgDy/XD9ZqT0amDRFC/tW6O79DMKatOVdSkhyGGYMM6on1Go1Dl+7gN2XjiLuThxSCi/oPDQp1MYwE55wN2+KNk6tEOIbiC5eLSA3khuocqLHY5gxzKgeS8m5h1//OY7YGydwKfs8slSXAXlB+Y5qBaxk3vCyao62zq3Qs1F7tHXx4a24qNZgmDHMiDTUajWO37iI6KsncPr2WVzPu4ACWbL2s9zKqKxgLfOGj3VztHMJQA+ftmjt2pABRwbBMGOYET1SfnEx9l05g4PXTuF8xnmkFl5GsfxmuVtxAfhvwHnB06oxWjv5oZt3W3T0bMpDlKR3DDOGGVGVZRXm448rcfgz+TQS7p7H7aLLKDa6XW70JABAbQpzeMLVvBH87Jujk0cAuvsGwMbMvOYLpzqLYcYwI6oWWYV52Hflb8TeOIuEjASkFFxBkdFNnYcohTCCicoVDqbe8LVpjDbOLfCcdwD8nNx5mJKeCMOMYUakN/nFxTiUGI/YG2dwPj0BN/MvI09c0z3IBABUFrCUecDF3AfN7ZugvZsfnvNpBSdL2xqtm6SHYcYwI6pRarUa59KSsT/xNP6+HY+k7Cu4W5qMUqO08ncw+S+Zyh5KI0+4WfqghUNTdHBviWcb+sHGnNfD0X0MM4YZUa2QVZiPQ4nncfTmefyTcQG38hORq74OGGfr7C+EEeQqRyiN3eFm4YVmDo3QzqU5nvX2g4MF/z3XNwwzhhlRrZZ4Nw0x187hdGoCLmdewu3CJBTgJmBUWOE0RipbWBm5w8W8IXxtGyHAqQk6efqhsYMzb7xcRzHMGGZEkiOEwD/pNxCbHI+zty/iStZVpBUkIx8pgDyn4glVFjCDGxxMPeBh3RDN7H3Q1rUJnvFoAqWZRc2tAFU7hhnDjKjOEELg6t07iE2Ox5nbl3A16ypSC64hV30TKqN7FZ6TE0IGI7UdrIxc4GTmAS8bLzR39EEb5yZo4+YDcxM+aaC2Y5gxzIjqhcyCXBy5/g9Op1zAxbuJuJmXjHvFt1Aku/3IQ5b3z805wEp+P+g8rb3QzMEbAc6+aOvmw5sy1xIMM4YZUb2mVqtxOeM2jt+4gHN3ruBqZhJSC24gpzQFJUZpum/l9V9CyGCksoOFUQPYK1zhZumORrYN4efkg7ZujeBp04Dn6GoIw4xhRkQVKFWpEH/nOk7evIyEjCtIyrqG2wXXkVN6G8Wy9EcGHQBArYCJcIS13BmOZq7wsPJAIzsvtHTyRmtXHzhaWtXMitQDDDOGGRE9AbVajct3U3Dq1hX8k56Eq5nJSM27iXslqSgUdwDjrMfPRGUFU9jDWu4ERzNnuFm5wcfOA80dG6K1izdcrR24Z1dJDDOGGRHpwb2CPMTduopzdxJx5W4ybuTeQHrhLWSXpqFEdgcwKn78TNSmMBb2sDRyhK2pM5wtXOFp7QpfO0/4OXnBz8kDFqam+l8ZCWCYMcyIqIYJIXAr5y7OpCbinzvJSMy8gZu5t5BReBs5pXdQLMsA5LmVmI8RZCprmMIOVsYOsFM0gLOFEzxsXOFj64ZmDTzh18ADlqZ1/6bODDOGGRHVQum5OTiXlox/7lxDYuZN3My9hbSCVGSXpKFAZPz3UgMdTynQRWUJE9jBUm4PGxNHOJo5wdXKGQ1tXOFr74bmjp7wVNpDLpfuTZ4ZZgwzIpKgUlUprt67jYS067h87yauZ6UgJS8VGYV3kFOSjkJxFyqjrMcPUvkvoTaBkdoGCpkSFnI72Jo6wMHcES6WDeBh4wwfO1c0tneFt50TTOQmel67qqvK97lxDdVERESPYSw3RlNHdzR1dK+wj0qlRuK9O7iYcRNX7t7E9exUpOSmIqMoDdnF6chX30UJMgF5PmRGJRBGGShEBgoB3C0GrhYDyAJw63/zvH+BuRWMoYS5kS2sje1hq3BAA4sGcLF0hKeNC7xsXdDI3gWu1spaubfHPTMiojoouygPlzNScDnjFq5l3cat7Nu4nZ+Ou4XpyCm9iwJVJkqQCSHPrfAuKroItTGM1NYwgQ3M5LawNraFUmEHR3NHOFs4wM2mAdytneCpdEJDpSOszUyfePQmDzMyzIiIKiW/uBhX797G1bupuJaVilu5aUjLu4OMwgxkFWcgX3UPRSITKll25UZrPkAIGTyMg7F75OdPVBsPMz4ttQo4twXw7ADYeRu6GiIivbEwNYW/iyf8XTwf2zerMBdX7qYi6d5t3Mi+g1s5d5CWn457hXeRVXIP+aWZKBJZKEUOIM+DTCZgJq+ZW4MxzHQ5/S3wy5T7/7+wEhdJEhHVA0ozK7Rza4x2bo0f27d4SSNkFd1F4eCQGqgMqH1n8WqDpMOGroCISNJM89PRQKWG582jNbI8hplOvNUMEVH1qJnvU4aZLrxvGhFR9aih71OGmU4MMyKi6sEwMxzumRERVQ/umdUSVw8augIiImm5fvyBNwwzA3pg4/880XBlEBFJ0Xcv/e//uWdmQA9ufFHJO1gTEdF9RdkPvKkjYRYZGQkfHx+YmZkhMDAQMTExj+x/8OBBBAYGwszMDL6+vvjyyy/1XaIOD2z8vDQDLJ+IiKpCr2G2adMmTJ06FXPnzsXp06fRtWtX9O7dG8nJyTr7JyYmok+fPujatStOnz6NOXPmYPLkydiyZYs+yyzvwT8kVMVAQWbNLp+ISKry72q/rwuHGZcuXYpx48Zh/PjxaNGiBZYtWwZPT0+sXLlSZ/8vv/wSDRs2xLJly9CiRQuMHz8eY8eOxccff6zPMh8vJc6wyycikooLux5qqJkw09u9GYuLi3Hy5EnMmjVLqz00NBSxsbE6pzly5AhCQ0O12sLCwrBmzRqUlJTAxKT8w+OKiopQVFSkeZ+dnV2uT5XJHsr49QOA+feA4lwgbiNwfiuQceX+XpuqBECdevAAkURV8KWpc89AR9tT9auob2X7VWWe1bw+TzVPGWDnBTTqAbR8EbDzAX7+VyXmXf30Fmbp6elQqVRwdnbWand2dkZqaqrOaVJTU3X2Ly0tRXp6OlxdXctNExERgXfffbf6Cgeg84e4yA4wMgHUlXvCKxFRvZBzC0g+Auz/oIIOEg+zMg8/lE0I8cgHtenqr6u9zOzZszF9+nTN++zsbHh6Pv5RBk9EXQI4NgU6vA54dgRMLQEjefk9OSKqeTofzVjBUZOn7lsbaqgFfdWlwO1zQPwO4Or+++8fVkP3oNBbmDk6OkIul5fbC0tLSyu391XGxcVFZ39jY2M4ODjonEahUEChUFRP0WV0BedLa+4/30zpyTuEEBGVcW8HtBt1f+BH4kFg85iHOkh8AIipqSkCAwMRHR2t1R4dHY3OnTvrnCYoKKhc/z179qB9+/Y6z5fpz0Mbv9scoNVgwLYhg4yISBcL+/vnzbrNMcji9Xp8bPr06fj666+xdu1aJCQkYNq0aUhOTsaECRMA3D9EOGrUKE3/CRMm4Nq1a5g+fToSEhKwdu1arFmzBjNmzNBnmeU9HFhdp+vuR0RE2p576Pu6hk7D6PWc2bBhw5CRkYFFixYhJSUF/v7+2LVrF7y8vAAAKSkpWtec+fj4YNeuXZg2bRpWrFgBNzc3LF++HC+99FJFi9CTh8JMXpN7hUREEmYk134v9dGMZd566y289dZbOj9bt25dubbg4GCcOnVKz1VVgYmFoSsgIpIWM1ugMPP+/9fQLQE5DE+XB/+SGLnVcHUQEUlRv0//9/9qhpkBPRBmPMRIRFQ1WjdrV9XIIhlmumgd4+XoRSKiqnnge1PNMKsdmGVERFXz4AhGnjMzoAf/knBobLg6iIikqHHI//6fhxkNqGzjtx8HmCkNWwsRkdSYWgId719PzMOMhlS28W3cDFsHEZFUyU3v/5eHGQ2o7IaavIEwEdGTKbt4mntmBlR2mPHhK9mJiKhyZP/9/uQ5MwMq+0tCxjAjInoiZTsDPMxoQNwzIyJ6OmWnaXiY0YC4Z0ZE9HRq+DCj3m80LEku/kBRNqD0MHQlRETS1LAj8Ow0wK1djSxOJkRVngle+2VnZ0OpVCIrKws2NjaGLoeIiJ5QVb7PeZiRiIgkj2FGRESSxzAjIiLJY5gREZHkMcyIiEjyGGZERCR5DDMiIpI8hhkREUkew4yIiCSPYUZERJLHMCMiIsljmBERkeQxzIiISPIYZkREJHkMMyIikjyGGRERSR7DjIiIJI9hRkREkscwIyIiyWOYERGR5DHMiIhI8hhmREQkeQwzIiKSPIYZERFJHsOMiIgkT69hdu/ePYSHh0OpVEKpVCI8PByZmZkV9i8pKcE777yDVq1awdLSEm5ubhg1ahRu3bqlzzKJiEji9Bpmw4cPR1xcHHbv3o3du3cjLi4O4eHhFfbPz8/HqVOn8H//9384deoUtm7diosXL+KFF17QZ5lERCRxMiGE0MeMExIS4Ofnh6NHj6Jjx44AgKNHjyIoKAj//PMPmjVrVqn5/PXXX+jQoQOuXbuGhg0blvu8qKgIRUVFmvfZ2dnw9PREVlYWbGxsqmdliIioxmVnZ0OpVFbq+1xve2ZHjhyBUqnUBBkAdOrUCUqlErGxsZWeT1ZWFmQyGWxtbXV+HhERoTmMqVQq4enp+bSlExGRxOgtzFJTU+Hk5FSu3cnJCampqZWaR2FhIWbNmoXhw4dXmMqzZ89GVlaW5nX9+vWnqpuIiKSnymG2cOFCyGSyR75OnDgBAJDJZOWmF0LobH9YSUkJXn75ZajVakRGRlbYT6FQwMbGRutFRET1i3FVJ5g4cSJefvnlR/bx9vbGmTNncPv27XKf3blzB87Ozo+cvqSkBEOHDkViYiL27dvHgCIiokeqcpg5OjrC0dHxsf2CgoKQlZWF48ePo0OHDgCAY8eOISsrC507d65wurIgu3TpEvbv3w8HB4eqlkhERPWM3s6ZtWjRAs8//zxee+01HD16FEePHsVrr72Gfv36aY1kbN68ObZt2wYAKC0txeDBg3HixAls2LABKpUKqampSE1NRXFxsb5KJSIiidPrdWYbNmxAq1atEBoaitDQUAQEBODbb7/V6nPhwgVkZWUBAG7cuIEdO3bgxo0baNOmDVxdXTWvqoyAJCKi+kVv15kZSlWuSyAiotqrVlxnRkREVFMYZkREJHkMMyIikjyGGRERSR7DjIiIJI9hRkREkscwIyIiyWOYERGR5DHMiIhI8hhmREQkeQwzIiKSPIYZERFJHsOMiIgkj2FGRESSxzAjIiLJY5gREZHkMcyIiEjyGGZERCR5DDMiIpI8hhkREUkew4yIiCSPYUZERJLHMCMiIsljmBERkeQxzIiISPIYZkREJHkMMyIikjyGGRERSR7DjIiIJI9hRkREkscwIyIiyWOYERGR5DHMiIhI8hhmREQkeQwzIiKSPIYZERFJHsOMiIgkj2FGRESSp9cwu3fvHsLDw6FUKqFUKhEeHo7MzMxKT//GG29AJpNh2bJlequRiIikT69hNnz4cMTFxWH37t3YvXs34uLiEB4eXqlpt2/fjmPHjsHNzU2fJRIRUR1grK8ZJyQkYPfu3Th69Cg6duwIAFi9ejWCgoJw4cIFNGvWrMJpb968iYkTJ+L3339H37599VUiERHVEXrbMzty5AiUSqUmyACgU6dOUCqViI2NrXA6tVqN8PBwzJw5Ey1btnzscoqKipCdna31IiKi+kVvYZaamgonJ6dy7U5OTkhNTa1wusWLF8PY2BiTJ0+u1HIiIiI05+SUSiU8PT2fuGYiIpKmKofZwoULIZPJHvk6ceIEAEAmk5WbXgihsx0ATp48ic8++wzr1q2rsM/DZs+ejaysLM3r+vXrVV0lIiKSuCqfM5s4cSJefvnlR/bx9vbGmTNncPv27XKf3blzB87Ozjqni4mJQVpaGho2bKhpU6lU+Pe//41ly5YhKSmp3DQKhQIKhaJqK0FERHVKlcPM0dERjo6Oj+0XFBSErKwsHD9+HB06dAAAHDt2DFlZWejcubPOacLDw9GzZ0+ttrCwMISHh+PVV1+taqlERFRP6G00Y4sWLfD888/jtddew1dffQUAeP3119GvXz+tkYzNmzdHREQEXnzxRTg4OMDBwUFrPiYmJnBxcXnk6EciIqrf9Hqd2YYNG9CqVSuEhoYiNDQUAQEB+Pbbb7X6XLhwAVlZWfosg4iI6jiZEEIYuojqlJ2dDaVSiaysLNjY2Bi6HCIiekJV+T7nvRmJiEjyGGZERCR5DDMiIpI8hhkREUkew4yIiCSPYUZERJLHMCMiIsljmBERkeQxzIiISPIYZkREJHkMMyIikjyGGRERSR7DjIiIJI9hRkREkscwIyIiyWOYERGR5BkbuoDqVvas0ezsbANXQkRET6Pse7wyz5Cuc2GWk5MDAPD09DRwJUREVB1ycnKgVCof2UcmKhN5EqJWq3Hr1i1YW1tDJpM90Tyys7Ph6emJ69evP/ZR3fUBt4c2bo//4bbQxu2h7Wm3hxACOTk5cHNzg5HRo8+K1bk9MyMjI3h4eFTLvGxsbPgL+QBuD23cHv/DbaGN20Pb02yPx+2RleEAECIikjyGGRERSR7DTAeFQoEFCxZAoVAYupRagdtDG7fH/3BbaOP20FaT26PODQAhIqL6h3tmREQkeQwzIiKSPIYZERFJHsOMiIgkj2FGRESSxzDTITIyEj4+PjAzM0NgYCBiYmIMXZLeRURE4JlnnoG1tTWcnJwwcOBAXLhwQauPEAILFy6Em5sbzM3N0a1bN5w/f95AFdeciIgIyGQyTJ06VdNW37bFzZs3MXLkSDg4OMDCwgJt2rTByZMnNZ/Xp+1RWlqKefPmwcfHB+bm5vD19cWiRYugVqs1fery9jh06BD69+8PNzc3yGQybN++Xevzyqx7UVERJk2aBEdHR1haWuKFF17AjRs3nq4wQVp++OEHYWJiIlavXi3i4+PFlClThKWlpbh27ZqhS9OrsLAwERUVJc6dOyfi4uJE3759RcOGDUVubq6mz4cffiisra3Fli1bxNmzZ8WwYcOEq6uryM7ONmDl+nX8+HHh7e0tAgICxJQpUzTt9Wlb3L17V3h5eYkxY8aIY8eOicTERPHHH3+Iy5cva/rUp+3x/vvvCwcHB/Hrr7+KxMREsXnzZmFlZSWWLVum6VOXt8euXbvE3LlzxZYtWwQAsW3bNq3PK7PuEyZMEO7u7iI6OlqcOnVKdO/eXbRu3VqUlpY+cV0Ms4d06NBBTJgwQautefPmYtasWQaqyDDS0tIEAHHw4EEhhBBqtVq4uLiIDz/8UNOnsLBQKJVK8eWXXxqqTL3KyckRTZo0EdHR0SI4OFgTZvVtW7zzzjvi2WefrfDz+rY9+vbtK8aOHavVNmjQIDFy5EghRP3aHg+HWWXWPTMzU5iYmIgffvhB0+fmzZvCyMhI7N69+4lr4WHGBxQXF+PkyZMIDQ3Vag8NDUVsbKyBqjKMrKwsAIC9vT0AIDExEampqVrbRqFQIDg4uM5um3/961/o27cvevbsqdVe37bFjh070L59ewwZMgROTk5o27YtVq9erfm8vm2PZ599Fnv37sXFixcBAH///TcOHz6MPn36AKh/2+NBlVn3kydPoqSkRKuPm5sb/P39n2r71Lm75j+N9PR0qFQqODs7a7U7OzsjNTXVQFXVPCEEpk+fjmeffRb+/v4AoFl/Xdvm2rVrNV6jvv3www84deoU/vrrr3Kf1bdtcfXqVaxcuRLTp0/HnDlzcPz4cUyePBkKhQKjRo2qd9vjnXfeQVZWFpo3bw65XA6VSoUPPvgAr7zyCoD69/vxoMqse2pqKkxNTWFnZ1euz9N8zzLMdHj4OWhCiCd+NpoUTZw4EWfOnMHhw4fLfVYfts3169cxZcoU7NmzB2ZmZhX2qw/bArj/jMD27dvjP//5DwCgbdu2OH/+PFauXIlRo0Zp+tWX7bFp0yZ899132LhxI1q2bIm4uDhMnToVbm5uGD16tKZffdkeujzJuj/t9uFhxgc4OjpCLpeX++sgLS2t3F8addWkSZOwY8cO7N+/X+u5cC4uLgBQL7bNyZMnkZaWhsDAQBgbG8PY2BgHDx7E8uXLYWxsrFnf+rAtAMDV1RV+fn5abS1atEBycjKA+vW7AQAzZ87ErFmz8PLLL6NVq1YIDw/HtGnTEBERAaD+bY8HVWbdXVxcUFxcjHv37lXY50kwzB5gamqKwMBAREdHa7VHR0ejc+fOBqqqZgghMHHiRGzduhX79u2Dj4+P1uc+Pj5wcXHR2jbFxcU4ePBgnds2ISEhOHv2LOLi4jSv9u3bY8SIEYiLi4Ovr2+92RYA0KVLl3KXaVy8eBFeXl4A6tfvBgDk5+eXe+qxXC7XDM2vb9vjQZVZ98DAQJiYmGj1SUlJwblz555u+zzx0JE6qmxo/po1a0R8fLyYOnWqsLS0FElJSYYuTa/efPNNoVQqxYEDB0RKSormlZ+fr+nz4YcfCqVSKbZu3SrOnj0rXnnllToz3PhxHhzNKET92hbHjx8XxsbG4oMPPhCXLl0SGzZsEBYWFuK7777T9KlP22P06NHC3d1dMzR/69atwtHRUbz99tuaPnV5e+Tk5IjTp0+L06dPCwBi6dKl4vTp05rLlyqz7hMmTBAeHh7ijz/+EKdOnRI9evTg0Hx9WLFihfDy8hKmpqaiXbt2muHpdRkAna+oqChNH7VaLRYsWCBcXFyEQqEQzz33nDh79qzhiq5BD4dZfdsWv/zyi/D39xcKhUI0b95crFq1Suvz+rQ9srOzxZQpU0TDhg2FmZmZ8PX1FXPnzhVFRUWaPnV5e+zfv1/nd8Xo0aOFEJVb94KCAjFx4kRhb28vzM3NRb9+/URycvJT1cXnmRERkeTxnBkREUkew4yIiCSPYUZERJLHMCMiIsljmBERkeQxzIiISPIYZkREJHkMMyIikjyGGRERSR7DjIiIJI9hRkREkvf/4b5/Dk9MxpAAAAAASUVORK5CYII=",
      "text/plain": [
       "<fipy.viewers.matplotlibViewer.matplotlib1DViewer.Matplotlib1DViewer at 0x130264910>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbMAAAGxCAYAAADoJViMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABahklEQVR4nO3deVwU9f8H8NeywHIvl9zI4Y2IB6aiGSoKeaWZR6moqZX19f5qef3UrL6klZkllqaYpWXmkaWZ5InhkQd5QJ4gHiCCct+7n98fftmvK4uCsiwDr+fjsY/az35m5j0D7ouZ+cyMTAghQEREJGFGhi6AiIjoaTHMiIhI8hhmREQkeQwzIiKSPIYZERFJHsOMiIgkj2FGRESSxzAjIiLJY5gREZHkMcyo1srJycHbb7+N0NBQNGjQADKZDAsXLtTZt1u3bpDJZJDJZDAyMoK1tTUaN26MIUOG4KeffoJara7UMrt16wZ/f/9qXAvd1q1bB5lMhqSkJL0vqzo96mfwOJGRkVi3bl259qSkJMhkMp2fEVUWw4xqrYyMDKxatQpFRUUYOHDgY/v7+vriyJEjiI2Nxfbt2zFr1iwUFBRgyJAh6NatG7KysvRfNFWoojBzdXXFkSNH0Ldv35oviuoMY0MXQFQRLy8v3Lt3DzKZDOnp6fj6668f2d/c3BydOnXSahs/fjyioqIwduxYvP7669i0aZM+S652BQUFMDc3N3QZeqVQKMr93IiqintmVGuVHTZ8Wq+++ir69OmDzZs349q1a5Wa5q+//kLXrl1hYWEBX19ffPjhh5pDlbm5ubC1tcUbb7xRbrqkpCTI5XJ89NFHmrajR4+iS5cuMDMzg5ubG2bPno2SkpJy03p7e6Nfv37YunUr2rZtCzMzM7z77rsAgHPnzmHAgAGws7ODmZkZ2rRpg2+++aZS67J582Z07NgRSqVSsz5jx47V6pOcnIyRI0fCyckJCoUCLVq0wCeffPLYw7MLFy7U+TN6+DCqt7c3zp8/j4MHD2p+rt7e3pptpusw4+HDhxESEgJra2tYWFigc+fO2Llzp87l7N+/H2+++SYcHR3h4OCAQYMG4datW5XaPlQ3MMyoXnjhhRcghEBMTMxj+6ampmLEiBEYOXIkduzYgd69e2P27Nn47rvvAABWVlYYO3YsNmzYUO7QZWRkJExNTTVhER8fj5CQEGRmZmLdunX48ssvcfr0abz//vs6l33q1CnMnDkTkydPxu7du/HSSy/hwoUL6Ny5M86fP4/ly5dj69at8PPzw5gxY7BkyZJHrsuRI0cwbNgw+Pr64ocffsDOnTsxf/58lJaWavrcuXMHnTt3xp49e/Dee+9hx44d6NmzJ2bMmIGJEyc+dntVxrZt2+Dr64u2bdviyJEjOHLkCLZt21Zh/4MHD6JHjx7IysrCmjVr8P3338Pa2hr9+/fXuXc9fvx4mJiYYOPGjViyZAkOHDiAkSNHVkvtJBGCSALu3LkjAIgFCxbo/Dw4OFi0bNmywul/++03AUAsXrz4kcsJDg4WAMSxY8e02v38/ERYWJjm/ZUrV4SRkZH49NNPNW0FBQXCwcFBvPrqq5q2YcOGCXNzc5GamqppKy0tFc2bNxcARGJioqbdy8tLyOVyceHCBa1lv/zyy0KhUIjk5GSt9t69ewsLCwuRmZlZ4fp8/PHHAsAj+8yaNUvnOr/55ptCJpNp1fPwz2DBggVC19dIVFRUufVr2bKlCA4OLtc3MTFRABBRUVGatk6dOgknJyeRk5OjaSstLRX+/v7Cw8NDqNVqreW89dZbWvNcsmSJACBSUlIqXG+qW7hnRvWCqMJj+1xcXNChQwettoCAAK1DlL6+vujXrx8iIyM18964cSMyMjK09mb279+PkJAQODs7a9rkcjmGDRumc9kBAQFo2rSpVtu+ffsQEhICT09PrfYxY8YgPz8fR44cqXBdnnnmGQDA0KFD8eOPP+LmzZvl+uzbtw9+fn7l1nnMmDEQQmDfvn0Vzl8f8vLycOzYMQwePBhWVlaadrlcjvDwcNy4cQMXLlzQmuaFF17Qeh8QEAAAlT6sTNLHMKN6oexLzc3N7bF9HRwcyrUpFAoUFBRotU2ZMgWXLl1CdHQ0AGDFihUICgpCu3btNH0yMjLg4uJSbn662oD7I/selpGRobO9bF0yMjIqXJfnnnsO27dvR2lpKUaNGgUPDw/4+/vj+++/r5b568O9e/cghKhSTQ//zBQKBQCU+5lR3cUwo3phx44dkMlkeO6556ptnj169IC/vz+++OILxMbG4tSpU/jXv/6l1cfBwQGpqanlptXVBkDnYAoHBwekpKSUay8b4ODo6PjIOgcMGIC9e/ciKysLBw4cgIeHB4YPH67Zo3ua+ZuZmQEAioqKtNrT09MfWdOj2NnZwcjI6KnWmeofhhnVeVFRUfjtt9/wyiuvoGHDhtU678mTJ2Pnzp2YPXs2nJ2dMWTIEK3Pu3fvjr179+L27duaNpVKVaVLBEJCQrBv375yo/PWr18PCwuLSg9rVygUCA4OxuLFiwEAp0+f1sw/Pj4ep06dKjd/mUyG7t27VzjPshGJZ86c0Wr/5ZdfdC6/MntKlpaW6NixI7Zu3arVX61W47vvvoOHh0e5Q7FEvM6MarXffvsNeXl5yMnJAXB/dOBPP/0EAOjTpw8sLCw0fQsKCnD06FHN/1+9ehXbt2/Hr7/+iuDgYHz55ZfVXt/IkSMxe/ZsHDp0CPPmzYOpqanW5/PmzcOOHTvQo0cPzJ8/HxYWFlixYgXy8vIqvYwFCxbg119/Rffu3TF//nzY29tjw4YN2LlzJ5YsWQKlUlnhtPPnz8eNGzcQEhICDw8PZGZm4rPPPoOJiQmCg4MBANOmTcP69evRt29fLFq0CF5eXti5cyciIyPx5ptvPjI4+vTpA3t7e4wbNw6LFi2CsbEx1q1bh+vXr5fr26pVK/zwww/YtGkTfH19YWZmhlatWumcb0REBHr16oXu3btjxowZMDU1RWRkJM6dO4fvv/++Wi7ZoDrGsONPiB7Ny8tLAND5enCkXNkoxLKXpaWl8PX1FYMHDxabN28WKpWqUsuraFTk6NGjhZeXl85pxowZI4yNjcWNGzd0fv7nn3+KTp06CYVCIVxcXMTMmTPFqlWrdI5m7Nu3r855nD17VvTv318olUphamoqWrdurTX6ryK//vqr6N27t3B3dxempqbCyclJ9OnTR8TExGj1u3btmhg+fLhwcHAQJiYmolmzZuKjjz4qt92gY0Tp8ePHRefOnYWlpaVwd3cXCxYsEF9//XW59UtKShKhoaHC2tpaANBsT12jGYUQIiYmRvTo0UNYWloKc3Nz0alTJ/HLL79o9SkbzfjXX39pte/fv18AEPv373/sNqK6QSZEFYZ5EZGW4uJieHt749lnn8WPP/5o6HKI6i0eZiR6Anfu3MGFCxcQFRWF27dvY9asWYYuiaheY5gRPYGdO3fi1VdfhaurKyIjI7WG4xNRzeNhRiIikjy9Ds0/dOgQ+vfvDzc3N8hkMmzfvv2x0xw8eBCBgYEwMzODr6+vXkagERFR3aLXMMvLy0Pr1q3xxRdfVKp/YmIi+vTpg65du+L06dOYM2cOJk+ejC1btuizTCIikrgaO8wok8mwbdu2Rz5k8Z133sGOHTuQkJCgaZswYQL+/vvvR95/joiI6rdaNQDkyJEjCA0N1WoLCwvDmjVrUFJSAhMTk3LTFBUVad1KR61W4+7du3BwcOCFlUREEiaEQE5ODtzc3GBk9OgDibUqzFJTU7XuLg4Azs7OKC0tRXp6us4bj0ZERGgeYEhERHXP9evX4eHh8cg+tSrMgPI3Wi07ClrRXtbs2bMxffp0zfusrCw0bNgQ169fh42Njf4KJSIivcrOzoanpyesra0f27dWhZmLi0u5u4mnpaXB2NhY52M5gPs3Ly173MODbGxsGGZERHVAZU4Z1aq75gcFBWmeDVVmz549aN++vc7zZURERICewyw3NxdxcXGIi4sDcH/ofVxcHJKTkwHcP0Q4atQoTf8JEybg2rVrmD59OhISErB27VqsWbMGM2bM0GeZREQkcXo9zHjixAmtZyGVndsaPXo01q1bh5SUFE2wAYCPjw927dqFadOmYcWKFXBzc8Py5cvx0ksv6bNMIiKSuDp3O6vs7GwolUpkZWXxnBnpnUqlQklJiaHLIJIsuVwOY2NjnefFqvJ9XqsGgBBJSW5uLm7cuIE69vcgUY2zsLCAq6truYfbVgXDjOgJqFQq3LhxAxYWFmjQoAEv0Cd6AkIIFBcX486dO0hMTESTJk0ee3F0RRhmRE+gpKQEQgg0aNAA5ubmhi6HSLLMzc1hYmKCa9euobi4GGZmZk80n1o1NJ9IarhHRvT0nnRvTGse1VAHERGRQTHMiIhI8hhmREQkeQwzIqJq9OKLL8LOzg6DBw/WtF2/fh3dunWDn58fAgICsHnzZgNWWL1ycnLwzDPPoE2bNmjVqhVWr15tkDo4mpGIqBpNnjwZY8eOxTfffKNpMzY2xrJly9CmTRukpaWhXbt26NOnDywtLQ1YafWwsLDAwYMHYWFhgfz8fPj7+2PQoEEV3hxeX7hnRkQG061bN0ydOlUv887IyICTkxOSkpL0Mv+KdO/evdwjS1xdXdGmTRsAgJOTE+zt7XH37t0K56HP7VLd5HI5LCwsAACFhYVQqVSaGwkMHjwYS5curZE6GGZE9cyYMWMgk8kgk8lgYmICZ2dn9OrVC2vXroVara7RWrZu3Yr33ntP8746v8QjIiLQv39/eHt7AwAOHTqE/v37w83NDTKZDNu3b9c5XWRkJHx8fGBmZobAwEDExMRUSz1lTpw4AbVaDU9Pz2qdryFlZmaidevW8PDwwNtvvw1HR0cAwPz58/HBBx8gOztb7zUwzIjqoeeffx4pKSlISkrCb7/9hu7du2PKlCno168fSktLa6wOe3v7Sj14saoKCgqwZs0ajB8/XtOWl5eH1q1b44svvqhwuk2bNmHq1KmYO3cuTp8+ja5du6J3795aN0QPDAyEv79/udetW7ceW1dGRgZGjRqFVatWPd0K1qDKrK+trS3+/vtvJCYmYuPGjbh9+zYAICAgAN7e3tiwYYP+CxV1TFZWlgAgsrKyDF0K1WEFBQUiPj5eFBQUGLqUKhs9erQYMGBAufa9e/cKAGL16tVCCCHUarVYvHix8PHxEWZmZiIgIEBs3rxZa5rg4GAxadIkMXPmTGFnZyecnZ3FggULtPps3rxZ+Pv7CzMzM2Fvby9CQkJEbm6uZvopU6Zo6gKg9Xr33XeFvb29KCws1JrnoEGDRHh4eIXruGXLFuHo6Fjh5wDEtm3byrV36NBBTJgwQautefPmYtasWRXOS5f9+/eLl156SautsLBQdO3aVaxfv/6x0z+4XYQQ4rfffhM2Njbim2++EUIIkZ2dLYYPHy4sLCyEi4uLWLp0ablpKprvxIkTxZQpU4Stra1wcnISX331lcjNzRVjxowRVlZWwtfXV+zatatK6/ugCRMmiB9//FHzfuHChaJr166PnKaif09V+T7nnhkRAQB69OiB1q1bY+vWrQCAefPmISoqCitXrsT58+cxbdo0jBw5EgcPHtSa7ptvvoGlpSWOHTuGJUuWYNGiRZqH7KakpOCVV17B2LFjkZCQgAMHDmDQoEE6b8782WefISgoCK+99hpSUlKQkpKCf//731CpVNixY4emX3p6On799Ve8+uqrFa7LoUOH0L59+yqtf3FxMU6ePInQ0FCt9tDQUMTGxlZpXg8TQmDMmDHo0aMHwsPDqzTtDz/8gKFDh2L9+vWa5z9Onz4df/75J3bs2IHo6GjExMTg1KlTlZrfN998A0dHRxw/fhyTJk3Cm2++iSFDhqBz5844deoUwsLCEB4ejvz8/ErN7/bt25rDiNnZ2Th06BCaNWum+bxDhw44fvw4ioqKqrTeVcXRjETVQAiBghKVQZZtbiKvtttqNW/eHGfOnEFeXh6WLl2Kffv2ISgoCADg6+uLw4cP46uvvkJwcLBmmoCAACxYsAAA0KRJE3zxxRfYu3cvevXqhZSUFJSWlmLQoEHw8vICALRq1UrnspVKJUxNTWFhYQEXFxdN+/DhwxEVFYUhQ4YAADZs2AAPDw9069atwvVISkqCm5tbldY9PT0dKpUKzs7OWu3Ozs5ITU2t9HzCwsJw6tQp5OXlwcPDA9u2bUNRURE2bdqEgIAAzbm6b7/9tsJtUSYyMhJz5szBzz//rHk2ZE5ODr755hts3LgRISEhAICoqKhKr2/r1q0xb948APcfkPzhhx/C0dERr732GoD757lWrlyJM2fOoFOnTo+d340bNzBu3DgIISCEwMSJExEQEKD53N3dHUVFRUhNTdX8DugDw4yoGhSUqOA3/3eDLDt+URgsTKvnn7IQAjKZDPHx8SgsLESvXr20Pi8uLkbbtm212h784gLuj9xLS0sDcP+LMyQkBK1atUJYWBhCQ0MxePBg2NnZVbqm1157Dc888wxu3rwJd3d3REVFaQaxVKSgoOCJb1j78HzLtkll/f677t+Dqg6u2bJlC27fvo3Dhw+jQ4cOmvarV6+ipKREq02pVGrtDT3Kgz8vuVwOBwcHrVAtC/Oyn+HjBAYGIi4ursLPy27EXdk9vSfFw4xEpJGQkAAfHx/NF+/OnTsRFxenecXHx+Onn37SmsbExETrvUwm00wvl8sRHR2N3377DX5+fvj888/RrFkzJCYmVrqmtm3bonXr1li/fj1OnTqFs2fPYsyYMY+cxtHREffu3av0Msqmkcvl5fbC0tLSyu2t1YQ2bdqgQYMGiIqK0josW/b/ukK3MnT9vB5sK5tvdY1sLbsEoUGDBtUyv4pwz4yoGpibyBG/KMxgy64O+/btw9mzZzFt2jT4+flBoVAgOTlZ65Dik5DJZOjSpQu6dOmC+fPnw8vLC9u2bcP06dPL9TU1NYVKVf5w7fjx4/Hpp5/i5s2b6Nmz52OHtbdt2xbfffddleo0NTVFYGAgoqOj8eKLL2rao6OjMWDAgCrNqzo0atQIn3zyCbp16wa5XK4ZhdmoUSOYmJjg+PHjmu2QnZ2NS5cuPfXPSh/OnTsHDw8PzXB9fWGYEVUDmUxWbYf6akLZOQyVSoXbt29j9+7diIiIQL9+/TBq1CjI5XLMmDED06ZNg1qtxrPPPovs7GzExsbCysoKo0ePrtRyjh07hr179yI0NBROTk44duwY7ty5gxYtWujs7+3tjWPHjiEpKQlWVlawt7eHkZERRowYgRkzZmD16tVYv379Y5cbFhaG2bNn4969e5pDmrm5ubh8+bKmT2JiIuLi4mBvb4+GDRsCuD+wIjw8HO3bt0dQUBBWrVqF5ORkTJgwoVLrW92aNm2K/fv3o1u3bpq7iFhbW2P06NGYOXMm7O3t4eTkhAULFsDIyKhWPpIoJiam3KAafZDOvz4iqja7d++Gq6srjI2NYWdnh9atW2P58uUYPXq05tlS7733HpycnBAREYGrV6/C1tYW7dq1w5w5cyq9HBsbGxw6dAjLli1DdnY2vLy88Mknn6B37946+8+YMQOjR4+Gn58fCgoKkJiYCG9vb9jY2OCll17Czp07MXDgwMcut1WrVmjfvj1+/PFHvPHGGwDuX6xcNogCgGbPcPTo0Vi3bh0AYNiwYcjIyMCiRYuQkpICf39/7Nq1S68DFx6nWbNm2Ldvn2YP7ZNPPsHSpUsxYcIE9OvXDzY2Nnj77bdx/fr1Jz5PqC+FhYXYtm1bhecRq5NMVPZAq0RkZ2dDqVQiKysLNjY2hi6H6qjCwkIkJiZq7hRB+terVy+0aNECy5cvr1T/Xbt2YcaMGTh37ly1PPyxNsvLy4O7uzs++eQTjBs3ztDlaKxYsQI///wz9uzZ88h+Ff17qsr3OffMiKhWu3v3Lvbs2YN9+/Y98u4dD+vTpw8uXbqEmzdv1qlbRwHA6dOn8c8//6BDhw7IysrCokWLAMAg5/YexcTEBJ9//nmNLIthRkS1Wrt27XDv3j0sXry40sPPy0yZMkVPVRnexx9/jAsXLmgGrsTExOh9kEVVvf766zW2LIYZEdVqNX3Xeylo27YtTp48aegyapW6fSCZiIjqBYYZERFJHsOMiIgkj2FGRESSxzAjIiLJY5gREZHkMcyIiEjyGGZERCR5DDMiIpI8hhkREUkew4yIiCSPYUZERJLHMCOiJ9KtWzdMnTq11s6vTEZGBpycnGrkhsX6WoeqGjx4MJYuXWroMmoUw4yoHunfvz969uyp87MjR45AJpPh1KlTNVzVfVu3bsV7770HoHpDISIiAv3794e3t7em7dChQ+jfvz/c3Nwgk8mwfft2ndNGRkZqHhhZ9pgVQxozZgxmzZr12H7z58/HBx98gOzs7BqoqnbQe5hV9Zdhw4YNaN26NSwsLODq6opXX30VGRkZ+i6TqF4YN24c9u3bh2vXrpX7bO3atWjTpg3atWtngMoAe3t7WFtbV+s8CwoKsGbNGowfP16rPS8vD61bt37kwz43bdqEqVOnYu7cuTh9+jS6du2K3r17Izk5uVprrCy1Wo2dO3dW6gGcAQEB8Pb2xoYNG2qgslpC6NEPP/wgTExMxOrVq0V8fLyYMmWKsLS0FNeuXdPZPyYmRhgZGYnPPvtMXL16VcTExIiWLVuKgQMHVnqZWVlZAoDIysqqrtUgKqegoEDEx8eLgoICQ5dSJSUlJcLZ2VksXLhQqz0vL09YW1uLzz//XAghhFqtFosXLxY+Pj7CzMxMBAQEiM2bN2tNExwcLKZMmSKEEKKwsFBMmjRJNGjQQCgUCtGlSxdx/Phxrf4qlUp8+OGHolGjRsLU1FR4enqK999/v9z8Ro8eLQBovRITE8U333wj7O3tRWFhodZ8Bw0aJMLDw3Wu75YtW4Sjo+MjtwkAsW3btnLtHTp0EBMmTNBqa968uZg1a1aF83pwmwghxG+//SZsbGzEN998o/l84sSJYsqUKcLW1lY4OTmJr776SuTm5ooxY8YIKysr4evrK3bt2lVu3ocOHRJOTk5CpVKJzZs3C39/f2FmZibs7e1FSEiIyM3N1eq/cOFC0bVr10eue21R0b+nqnyf6zXMqvrL8NFHHwlfX1+ttuXLlwsPD49KL5NhRjWh3D8+tVqIolzDvNTqKtU+c+ZM4e3tLdQPTLdu3TqhUCjE3bt3hRBCzJkzRzRv3lzs3r1bXLlyRURFRQmFQiEOHDigmebBL+7JkycLNzc3sWvXLnH+/HkxevRoYWdnJzIyMjT93377bWFnZyfWrVsnLl++LGJiYsTq1avLzS8zM1MEBQWJ1157TaSkpIiUlBRRWloq8vPzhVKpFD/++KNmmjt37ghTU1Oxb98+nes6ZcoU8fzzzz9ye+gKs6KiIiGXy8XWrVu12idPniyee+65Cuf14Db5/vvvhbW1tdi+fbvW59bW1uK9994TFy9eFO+9954wMjISvXv3FqtWrRIXL14Ub775pnBwcBB5eXla854xY4YYN26cuHXrljA2NhZLly4ViYmJ4syZM2LFihUiJydHq/+uXbuEQqEoF/61UXWEmd6eNF1cXIyTJ0+WO74bGhqK2NhYndN07twZc+fOxa5du9C7d2+kpaXhp59+Qt++fStcTlFREYqKijTv69MxYqpFSvKB/7gZZtlzbgGmlpXuPnbsWHz00Uc4cOAAunfvDuD+IcZBgwbBzs4OeXl5WLp0Kfbt24egoCAAgK+vLw4fPoyvvvoKwcHBWvPLy8vDypUrsW7dOvTu3RsAsHr1akRHR2PNmjWYOXMmcnJy8Nlnn+GLL77A6NGjAQCNGjXCs88+W64+pVIJU1NTWFhYwMXFRdNubm6O4cOHIyoqCkOGDAFw/7SEh4cHunXrpnNdk5KS4OZW9Z9Leno6VCoVnJ2dtdqdnZ2Rmpr62OkjIyMxZ84c/Pzzz5ptXKZ169aYN28eAGD27Nn48MMP4ejoiNdeew3A/fNdK1euxJkzZ9CpUyfNdDt27MDHH3+MlJQUlJaWYtCgQfDy8gIAtGrVqlwN7u7uKCoqQmpqqqZfXaa3MHuSX4bOnTtjw4YNGDZsGAoLC1FaWooXXngBn3/+eYXLiYiIwLvvvluttRPVZc2bN0fnzp2xdu1adO/eHVeuXEFMTAz27NkDAIiPj0dhYSF69eqlNV1xcTHatm1bbn5XrlxBSUkJunTpomkzMTFBhw4dkJCQAABISEhAUVERQkJCnqr21157Dc888wxu3rwJd3d3REVFYcyYMZDJZDr7FxQUwMzM7ImX9/B8hRAVLqvMli1bcPv2bRw+fBgdOnQo93lAQIDm/+VyORwcHLTCqOw7My0tTdOWkJCAGzduoGfPnjA1NUVISAhatWqFsLAwhIaGYvDgwbCzs9Najrm5OQAgPz+/kmsrbXoLszJV+WWIj4/H5MmTMX/+fISFhSElJQUzZ87EhAkTsGbNGp3TzJ49G9OnT9e8z87OhqenZ/WtAFFlmFjc30My1LKraNy4cZg4cSJWrFiBqKgoeHl5aYJGrVYDAHbu3Al3d3et6RQKRbl5CSEAPPrfetkX69Nq27YtWrdujfXr1yMsLAxnz57FL7/8UmF/R0dH3Lt3r8rLcXR0hFwuL/eHd1paWrk/0B/Wpk0bnDp1ClFRUXjmmWfKbRcTExOt9zKZTKutrH/ZzwG4v1fWq1cvzXaMjo5GbGws9uzZg88//xxz587FsWPH4OPjo5nm7t27AIAGDRpUdrUlTW+jGZ/klyEiIgJdunTBzJkzERAQgLCwMERGRmLt2rVISUnROY1CoYCNjY3Wi6jGyWT3D/UZ4vWYPQVdhg4dCrlcjo0bN+Kbb77Bq6++qvkS9fPzg0KhQHJyMho3bqz10vWHYuPGjWFqaorDhw9r2kpKSnDixAm0aNECANCkSROYm5tj7969larP1NQUKpVK52fjx49HVFQU1q5di549ez7yj9e2bdsiPj6+Ust8ePmBgYGIjo7Wao+Ojkbnzp0fOW2jRo2wf/9+/Pzzz5g0aVKVl63Lzz//jBdeeEHzXiaToUuXLnj33Xdx+vRpmJqaYtu2bVrTnDt3Dh4eHnB0dKyWGmo7ve2ZPfjL8OKLL2rao6OjKxxamp+fD2Nj7ZLkcjmA//31R0RPz8rKCsOGDcOcOXOQlZWFMWPGaD6ztrbGjBkzMG3aNKjVajz77LPIzs5GbGwsrKysNOe8ylhaWuLNN9/EzJkzYW9vj4YNG2LJkiXIz8/HuHHjAABmZmZ455138Pbbb8PU1BRdunTBnTt3cP78eU2fB3l7e+PYsWNISkqClZUV7O3tYWR0/2/vESNGYMaMGVi9ejXWr1//yPUMCwvD7Nmzce/ePa3DcLm5ubh8+bLmfWJiIuLi4jT1A8D06dMRHh6O9u3bIygoCKtWrUJycjImTJjw2O3btGlT7N+/H926dYOxsTGWLVv22GkqkpaWhr/++ktzLdyxY8ewd+9ehIaGwsnJCceOHcOdO3c0fziUiYmJQWho6BMvV3L0MDBFo2xo/po1a0R8fLyYOnWqsLS0FElJSUIIIWbNmqU1pDYqKkoYGxuLyMhIceXKFXH48GHRvn170aFDh0ovk6MZqSZIdWj+g2JjYwUAERoaWu4ztVotPvvsM9GsWTNhYmIiGjRoIMLCwsTBgwc1fR4cuVdQUCAmTZokHB0dHzk0//333xdeXl7CxMRENGzYUPznP//ROb8LFy6ITp06CXNzc83Q/AeFh4frHKavS6dOncSXX36p1bZ///5yw/8BiNGjR2v1W7FihfDy8hKmpqaiXbt2Wuuvy8ND8+Pj44WTk5OYPn26zs+FEMLLy0t8+umnWm14YITl119/Lbp06aI1z7CwMM1lEE2bNtVcUlGmoKBA2NjYiCNHjjyy3tqi1g/NF+LRvwyjR48WwcHBWv2XL18u/Pz8hLm5uXB1dRUjRowQN27cqPTyGGZUE+pCmElZz549xaRJkyrVd+fOnaJFixZCpVLpuSr96N+/v1i8eHGVpvniiy9Er1699FRR9auOMJMJUbeO32VnZ0OpVCIrK4vnz0hvCgsLkZiYqLm7DdWMu3fvYs+ePRgxYgTi4+PRrFmzSk332WefYdCgQZIcHLZkyRK88sorVap91apVCA4OrvT2MbSK/j1V5ftc76MZiYiqS7t27XDv3j0sXry4Sl/UU6ZM0WNV+vX2229XeZrXX39dD5XUbgwzIpKMmrjzPUkT75pPRESSxzAjIiLJY5gREZHkMcyIiEjyGGZERCR5DDMiIpI8hhkREUkew4yIiCSPYUZERJLHMCOiemXNmjV1+tEogwcPxtKlSw1dRo1jmBHVI2PGjIFMJtM83djX1xczZsxAXl6eoUurEUVFRZg/fz7+7//+T6t9y5YtmoeS+vn5lXvQpS5nz55FcHAwzM3N4e7ujkWLFmk9d/Hw4cPo0qULHBwcYG5ujubNm+PTTz995DyTkpI0Px+ZTAY7Ozs899xzOHjwYKXXcf78+fjggw+QnZ1d6WnqAoYZUT3z/PPPIyUlBVevXsX777+PyMhIzJgxQ2ffkpKSGq7u8Z6mpi1btsDKygpdu3bVtB05cgTDhg1DeHg4/v77b4SHh2Po0KE4duxYhfPJzs5Gr1694Obmhr/++guff/45Pv74Y609IktLS0ycOBGHDh1CQkIC5s2bh3nz5mHVqlWPrfOPP/5ASkoKDh48CBsbG/Tp0weJiYmVWseAgAB4e3tjw4YNlepfZ+jl4TQGxOeZUU14+PlLarVa5BXnGeSlVqsrXffo0aPFgAEDtNrGjx8vXFxchBBCLFiwQLRu3VqsWbNG+Pj4CJlMJtRqtcjMzBSvvfaaaNCggbC2thbdu3cXcXFxmnnExcWJbt26CSsrK2FtbS3atWsn/vrrLyGEEElJSaJfv37C1tZWWFhYCD8/P7Fz504hxP0H8iqVSq16tm3bJh78anrSmnTp37+/mDFjhlbb0KFDxfPPP6/VFhYWJl5++eUK5xMZGSmUSqXWw0EjIiKEm5vbI38eL774ohg5cmSFnycmJgoA4vTp05q2GzduCADiyy+/FOnp6eLll18W7u7uwtzcXPj7+4uNGzeWm8/ChQtF165dK1xObVMdzzPjXfOJqkFBaQE6buxokGUfG34MFiYWTzy9ubm51t7O5cuX8eOPP2LLli2Qy+UAgL59+8Le3h67du2CUqnEV199hZCQEFy8eBH29vYYMWIE2rZti5UrV0IulyMuLg4mJiYAgH/9618oLi7GoUOHYGlpifj4eFhZWVWpxiepSZeYmBiMGDFCq+3IkSOYNm2aVltYWBiWLVtWYT1HjhxBcHAwFAqF1jSzZ89GUlISfHx8yk1z+vRpxMbG4v3336/sagMALCzu/2xLSkpQWFiIwMBAvPPOO7CxscHOnTsRHh4OX19fdOz4v9+/Dh06ICIiAkVFRVo11mUMM6J67Pjx49i4cSNCQkI0bcXFxfj222/RoEEDAMC+fftw9uxZpKWlab4YP/74Y2zfvh0//fQTXn/9dSQnJ2PmzJlo3rw5AKBJkyaa+SUnJ+Oll15Cq1atAAC+vr5VrvNJanpYZmYmMjMz4ebmptWempoKZ2dnrTZnZ2ekpqZWWE9qaiq8vb3LTVP22YNh5uHhgTt37qC0tBQLFy7E+PHjK73eeXl5mD17NuRyOYKDg+Hu7q51SHjSpEnYvXs3Nm/erBVm7u7uKCoqQmpqKry8vCq9PCljmBFVA3NjcxwbXvE5Fn0vuyp+/fVXWFlZobS0FCUlJRgwYAA+//xzzedeXl6a0ACAkydPIjc3Fw4ODlrzKSgowJUrVwAA06dPx/jx4/Htt9+iZ8+eGDJkCBo1agQAmDx5Mt58803s2bMHPXv2xEsvvYSAgIAq1fwkNT2soKAAAHQ+GVwmk2m9F0KUa6vMNLraY2JikJubi6NHj2LWrFlo3LgxXnnllUfOu3PnzjAyMkJ+fj5cXV2xbt06tGrVCiqVCh9++CE2bdqEmzdvoqioCEVFRbC0tNSa3tz8/u9Efn7+I5dTlzDMiKqBTCZ7qkN9Nal79+5YuXIlTExM4ObmpjkcWObhL0a1Wg1XV1ccOHCg3LxsbW0BAAsXLsTw4cOxc+dO/Pbbb1iwYAF++OEHvPjiixg/fjzCwsKwc+dO7NmzBxEREfjkk08wadIkGBkZaY0ABHQP8HiSmh7m4OAAmUyGe/fuabW7uLiU2wtLS0srt7dWmWkAlJuubC+tVatWuH37NhYuXPjYMNu0aRP8/Pxga2urFdiffPIJPv30UyxbtgytWrWCpaUlpk6diuLiYq3p7969CwBafwDUdRzNSFTPWFpaonHjxvDy8ioXZLq0a9cOqampMDY2RuPGjbVejo6Omn5NmzbFtGnTsGfPHgwaNAhRUVGazzw9PTFhwgRs3boV//73v7F69WoA979sc3JytC4NiIuLq7aaHmRqago/Pz/Ex8drtQcFBSE6Olqrbc+ePejcuXOFyw8KCsKhQ4e0QmTPnj1wc3Mrd/jxQUIIFBUVPXb9PD090ahRo3J7njExMRgwYABGjhyJ1q1bw9fXF5cuXSo3/blz5+Dh4VHhtqiLGGZE9Eg9e/ZEUFAQBg4ciN9//x1JSUmIjY3FvHnzcOLECRQUFGDixIk4cOAArl27hj///BN//fUXWrRoAQCYOnUqfv/9dyQmJuLUqVPYt2+f5rOOHTvCwsICc+bMweXLl7Fx40asW7fuqWuqSFhYGA4fPqzVNmXKFOzZsweLFy/GP//8g8WLF+OPP/7A1KlTNX2++OILrfOKw4cPh0KhwJgxY3Du3Dls27YN//nPfzB9+nTNYcYVK1bgl19+waVLl3Dp0iVERUXh448/xsiRIyu76ctp3LgxoqOjERsbi4SEBLzxxhs6z+3FxMTU6QvDddLHMEtD4tB8qgkVDSWu7XQNzX9Q2TD4h2VnZ4tJkyYJNzc3YWJiIjw9PcWIESNEcnKyKCoqEi+//LLw9PQUpqamws3NTUycOFGzbSZOnCgaNWokFAqFaNCggQgPDxfp6emaeW/btk00btxYmJmZiX79+olVq1bpHJpflZoqkpCQIMzNzUVmZqZW++bNm0WzZs2EiYmJaN68udiyZUu57eLl5aXVdubMGdG1a1ehUCiEi4uLWLhwodaw/OXLl4uWLVsKCwsLYWNjI9q2bSsiIyOFSqWqsD5dQ/MflJGRIQYMGCCsrKyEk5OTmDdvnhg1apTWz7SgoEDY2NiII0eOVLic2qY6hubLhHjogLXEZWdnQ6lUIisrCzY2NoYuh+qowsJCJCYmwsfHR+eAAqq9hg4dirZt22L27NmGLkUvVqxYgZ9//hl79uwxdCmVVtG/p6p8n/MwIxHVKx999FGVr3OTEhMTE63RqfUFRzMSUb3i5eWFSZMmGboMvdF1jV19wD0zIiKSPIYZ0VOoY6eciQyiOv4dMcyInkDZ/QEfvliViKqu7E4llbnusSI8Z0b0BIyNjWFhYYE7d+7AxMQERkb8u5CoqoQQyM/PR1paGmxtbTV/JD4JhhnRE5DJZHB1dUViYiKuXbtm6HKIJM3W1hYuLi5PNQ+GGdETMjU1RZMmTXiokegpmJiYPNUeWRmGGdFTMDIy4kXTRLUAD/QTEZHkMcyIiEjyGGZERCR5DDMiIpI8hhkREUkew4yIiCRP72EWGRmpeUZNYGAgYmJiHtm/qKgIc+fOhZeXFxQKBRo1aoS1a9fqu0wiIpIwvV5ntmnTJkydOhWRkZHo0qULvvrqK/Tu3Rvx8fFo2LChzmmGDh2K27dvY82aNWjcuDHS0tJQWlqqzzKJiEji9Pqk6Y4dO6Jdu3ZYuXKlpq1FixYYOHAgIiIiyvXfvXs3Xn75ZVy9ehX29vZPtEw+aZqIqG6oFU+aLi4uxsmTJxEaGqrVHhoaitjYWJ3T7NixA+3bt8eSJUvg7u6Opk2bYsaMGSgoKKhwOUVFRcjOztZ6ERFR/aK3w4zp6elQqVRwdnbWand2dkZqaqrOaa5evYrDhw/DzMwM27ZtQ3p6Ot566y3cvXu3wvNmERERePfdd6u9fiIikg69DwCRyWRa74UQ5drKqNVqyGQybNiwAR06dECfPn2wdOlSrFu3rsK9s9mzZyMrK0vzun79erWvAxER1W562zNzdHSEXC4vtxeWlpZWbm+tjKurK9zd3aFUKjVtLVq0gBACN27cQJMmTcpNo1AooFAoqrd4IiKSFL3tmZmamiIwMBDR0dFa7dHR0ejcubPOabp06YJbt24hNzdX03bx4kUYGRnBw8NDX6USEZHE6fUw4/Tp0/H1119j7dq1SEhIwLRp05CcnIwJEyYAuH+IcNSoUZr+w4cPh4ODA1599VXEx8fj0KFDmDlzJsaOHQtzc3N9lkpERBKm1+vMhg0bhoyMDCxatAgpKSnw9/fHrl274OXlBQBISUlBcnKypr+VlRWio6MxadIktG/fHg4ODhg6dCjef/99fZZJREQSp9frzAyB15kREdUNteI6MyIioprCMCMiIsljmBERkeQxzIiISPIYZkREJHkMMyIikjyGGRERSR7DjIiIJI9hRkREkscwIyIiyWOYERGR5DHMiIhI8hhmREQkeQwzIiKSPIYZERFJHsOMiIgkj2FGRESSxzAjIiLJY5gREZHkMcyIiEjyGGZERCR5DDMiIpI8hhkREUkew4yIiCSPYUZERJLHMCMiIsljmBERkeQxzIiISPIYZkREJHkMMyIikjyGGRERSR7DjIiIJI9hRkREkscwIyIiyWOY6XAn7x6Wn1qOEnWJoUshIqJKYJg9RC3UGLh1JFafXY1nvx2Ao8mXDF0SERE9BsPsIUYyI9gV9wEA5OM6xv8xHG/98hkKSooNXBkREVWEYabDT6MmYKr/RzBTN4RMXoiYu1+j8/oXsCHugKFLIyIiHWRCCGHoIqpTdnY2lEolsrKyYGNj81TzKlWVYu7er7HrxlpAXgAAcDHqgs/C5sHPyaM6yiUiogpU5ftc73tmkZGR8PHxgZmZGQIDAxETE1Op6f78808YGxujTZs2+i3wEYzlxlgcOgE7XvwFDU26QwgZUtV/YujOgXhjxyfIKy40WG1ERPQ/eg2zTZs2YerUqZg7dy5Onz6Nrl27onfv3khOTn7kdFlZWRg1ahRCQkL0WV6l+dg5Y+fw5Xivw1dQqLwhMypC7L116PJtP6w8tsvQ5RER1Xt6PczYsWNHtGvXDitXrtS0tWjRAgMHDkRERESF07388sto0qQJ5HI5tm/fjri4uEovszoPM+pSqlJhwf512JH8NSDPBQDYifb4qMdcdGzYuNqXR0RUX9WKw4zFxcU4efIkQkNDtdpDQ0MRGxtb4XRRUVG4cuUKFixYUKnlFBUVITs7W+ulT8ZyOT7oOQ67B/+KpubPQwgZ7slOYNzeIXhl80Kk5ep3+UREVJ7ewiw9PR0qlQrOzs5a7c7OzkhNTdU5zaVLlzBr1ixs2LABxsbGlVpOREQElEql5uXp6fnUtVeGu40Dtgz9CJ91XQ8r0Qwyo1Kcy9+CkB/7YMG+dVCpVTVSBxER1cAAEJlMpvVeCFGuDQBUKhWGDx+Od999F02bNq30/GfPno2srCzN6/r1609dc1WENGqDP0f9iLGNF8JI5QDIs7D1+ifouG4ANp87VKO1EBHVV3o7Z1ZcXAwLCwts3rwZL774oqZ9ypQpiIuLw8GDB7X6Z2Zmws7ODnK5XNOmVqshhIBcLseePXvQo0ePxy5X3+fMHiWnqADTdn+Bo3c3QWZUBABwMuqEj3rMRjt33xqthYhI6mrFOTNTU1MEBgYiOjpaqz06OhqdO3cu19/GxgZnz55FXFyc5jVhwgQ0a9YMcXFx6Nixo75KrTbWCnN8PWAmNvXZBnd5NwghQ5r6KEbteQmvbF6I1JwsQ5dIRFQnVe7E1BOaPn06wsPD0b59ewQFBWHVqlVITk7GhAkTANw/RHjz5k2sX78eRkZG8Pf315reyckJZmZm5dpru5bOntg98nP8+s8JvHckAvlGF3Eufwt6bf4dIS4j8Z9e42FhojB0mUREdYZez5kNGzYMy5Ytw6JFi9CmTRscOnQIu3btgpeXFwAgJSXlsdecSVm/5u1xZPRmvNp4PuSqBoA8F3vvfImg9b2xJOZHqNVqQ5dIRFQn8HZWNaSgpAhz//ga0Snfaa5PM1P7YFrgNAwP6G7g6oiIap+qfJ8zzGpYak4mZvy+HHE5P0NmdP9O/LZojfldZqBX4zaGLY6IqBZhmNXiMCtz/vZ1vP3HUlwr2QeZTA0hZHCTP4sPuk3HM568kwgREcNMAmFW5lBiPBbEfIR0cQIAIIQcPqY98GHINLR0rpkLwImIaiOGmYTCrMz2+FgsOf4ZcmTxAAChNkYLy+expOcU+Ni7GLg6IqKaxzCTYJiV+ebUH1gRtwIF8sv3G9SmCLDphw9DJsHT1tGwxRER1SCGmYTDDLh/55PIY7uwNn4lSoz/e+mC2gxtbF7A4p4T4aa0M2yBREQ1gGEm8TAro1Kp8Wnsdmy4uAqlxjf/22iBtsr++CDkTXjaOhi2QCIiPWKY1ZEwK1OqUuHjPzdj0+WvUSq/DQAQajO0tumL//R4C152TgaukIio+jHM6liYlSkpLcXHsT/ix8vrUCpPAQAItSn8rXvj/e7/QmMHVwNXSERUfRhmdTTMypSqVPg0diu+vxiFEuP7j7wRamM0Me+Jd4MnIsDVy8AVEhE9PYZZHQ+zMiqVGsuP7sB3/6xBsXESAECo5Who2g1zn30TXbyaGbZAIqKnwDCrJ2FWRq1WY/WJ37H2/GrkG10CAAghQwOjZzC1/esY4Ff7H59DRPQwhlk9C7MHbfx7PyLjViMLZzVtVuoWGOs/BuMCn4eRkd4fLk5EVC0YZvU4zMr8fvEUPjm2CrdURyCT3X/UjKnKAwN9R2BGlyEw5/PUiKiWY5gxzDRO3LiC9w9/hcsFezV36ZeV2qOr0yD8X7cxcLFWGrhCIiLdGGYMs3IS76bh3YNrcPLeDs3z1ITKHM0te2J2l/EI9PA1cIVERNoYZgyzCmUW5OE/Md9hz41NUMnvAACEMEIDo0C80WY0hrV6DjKZzMBVEhExzBhmlVCqKsWXf+3EhoQNyDVK0LQrVA3R32cYZnQeBksFz6sRkeEwzBhmVfL7pdNYdmwtrpcchsyo9H6jyhptbftgTtexaN7AzbAFElG9xDBjmD2RS+mp+E9MFE7c2wnIswDcvwjbWd4R41uPwLBWz3JoPxHVGIYZw+yp5JcUYemfP+HnxB9QaJSkaTcudUdPjxfxzrOvwNGS25aI9IthxjCrNjsSjuHLU98iufjP/x2CVCvQyDwYk54ZhZBGrQ1bIBHVWQwzhlm1u56ZgQ8Pf4vDt3+B2jhN026uaoK+3oMwrfNLsDEzN2CFRFTXMMwYZnqjUqmx+sQebPxnE+6KU5q7i0BlhSYWwZj4zAj0aNTKsEUSUZ3AMGOY1Yhzqcn4+Mh6nLq3G+K/A0YAwEzli16e/fHvzsPgYGltwAqJSMoYZgyzGlVUWoLVJ3Zh84UtyBB//29vTa2Al6IzXm09FC+26MyRkERUJQwzhpnBJKTdwNIjG3E8fTfUxnc07cYqF3R0fB7Tgl5BM163RkSVwDBjmBmcSqXG+tP7sTFhM1JKj0NmVALg/q2z7BCAPj798FbH/lCaWRi4UiKqrRhmDLNa5Vb2XXx6ZDP23/wVRfKk/32gMoeXWRBG+g/CUP+uPAxJRFoYZgyzWuvA1TP46tRmnM8+ACHP1LQbqRzQxjYEbwQORmevFoYrkIhqDYYZw6zWKyktxfq4vdiUsB23So9rnrUGAAqVD7q4hGFih5fQxNHFgFUSkSExzBhmkpKWm43Pj27H3uu/IVt2HjLZ/V9JIYxgI/zQwzMUb3UYADcbewNXSkQ1iWHGMJOshLQbWHF8M46kRaNYfl3TLtRyOBgFoLfP83jjmRdgZ25lwCqJqCYwzBhmdcKf1+Kx+tQ2xN09AJVxqqZdqE3gbNwWz3v3xuvP9OGISKI6imHGMKtT1Go19lz6G+vObEN8dgyEcbrmM6FWwEUeiOd9wjC+/fOwNWewEdUVDDOGWZ2lUqmxPeEYvj//Cy7mxWiNiBRqBRrIW6OXVwjGB/aFk5XScIUS0VNjmDHM6gWVWoUt8bHYFP8LLuUeeSjYjGFv1BLdPEIwrl1feNk5Ga5QInoiDDOGWb2jVqux45/j+OH8TiRk/6l1Ky0hjGAlmqKjUzBebdsPbdy8DVcoEVVaVb7P9X7LhcjISPj4+MDMzAyBgYGIiYmpsO/WrVvRq1cvNGjQADY2NggKCsLvv/+u7xKpDjAyMsJAv074Ych7OP3qH/goaD3aWA2DcakHZDI18oz+wb70rxAe3R/t1vTDqK0f4Nd/TkCtVhu6dCKqBnrdM9u0aRPCw8MRGRmJLl264KuvvsLXX3+N+Ph4NGzYsFz/qVOnws3NDd27d4etrS2ioqLw8ccf49ixY2jbtm2llsk9M3rYseSL+ObvX3Ei/RAKjK5ofSZT2cLbvD16N+qBEQEhsOHISKJao9YcZuzYsSPatWuHlStXatpatGiBgQMHIiIiolLzaNmyJYYNG4b58+dXqj/DjB7lUvpNRJ3+DYdvHcRd9XnNDZABAGoT2Bv5o7NrV4S36Q0/Jw/DFUpEVfo+N9ZXEcXFxTh58iRmzZql1R4aGorY2NhKzUOtViMnJwf29hXf+aGoqAhFRUWa99nZ2U9WMNULTRzd8Z9e4wGMR2ZBHr6N24vdV/chufAkYJyJuziNX1NO49eU5TApbQg/2054oUkPDGjRAQoTE0OXT0QV0FuYpaenQ6VSwdnZWavd2dkZqampFUyl7ZNPPkFeXh6GDh1aYZ+IiAi8++67T1Ur1U+25paYFPQCJgW9AJVKjZ0XTmLLP9E4n3kERfIklBgn4+/cZPx9+ke8d8ICTiYB6OzWGSMDQtGsgbuhyyeiB+gtzMrIZDKt90KIcm26fP/991i4cCF+/vlnODlVPKx69uzZmD59uuZ9dnY2PD09n7xgqpfkciO84PcMXvB7BgBwMf0mvv17N/68eRh3VOcAeT7S1Eex/cZRbL+xFKYqdzS1aY9Q32AMadkVVgozA68BUf2mtzBzdHSEXC4vtxeWlpZWbm/tYZs2bcK4ceOwefNm9OzZ85F9FQoFFArFU9dL9KCmju54L2QcgHEoKCnC9vgj+PXSAfyT9ReK5NdRLL+Jc3k3ce7sz/jkb1PYylqgbYOOeKlFCJ7zbs5nsxHVML0PAAkMDERkZKSmzc/PDwMGDKhwAMj333+PsWPH4vvvv8fAgQOrvEwOACF9u3AnBd/9HY2jKbFILTkDyHO0PpeV2sFV4Y9Orp0wzL87/Jx5pIDoSdSa0YxlQ/O//PJLBAUFYdWqVVi9ejXOnz8PLy8vzJ49Gzdv3sT69esB3A+yUaNG4bPPPsOgQYM08zE3N4dSWblbEzHMqCaVqlTYc/k0fr6wH2cyjiFHdhkymUqrj7HKGV4WrfGcRxCG+HeDp62jgaolkpZaE2bA/YumlyxZgpSUFPj7++PTTz/Fc889BwAYM2YMkpKScODAAQBAt27dcPDgwXLzGD16NNatW1ep5THMyJAyC3Kx+VwM9ib9iUs5cSgyStY8nw0AhJBBofaAr2UbPNcwCEP8g+Fizd9TIl1qVZjVNIYZ1SbX7qXjh7MHcPjGEVwvOKP1KBvg/q22zFRe8LFuha6eHfBSy65w50NIiQAwzBhmVGvFp93AT+cP4GjKMdwsPAO1/K7W5/f33NzhZemPILdn8KJfVzR2cDVQtUSGxTBjmJFExKVcwfaEwzieegK3CuKhMk4r10de6gQ3s5Zo59wO/Zs+iw6ejSp1eQuR1DHMGGYkUfG3b2BbQgyO3jqBGwXnUSK/pXXODQBQagsH4yZo6RCA7l7PoHfT9rDk5SlUBzHMGGZUR1zPTMe2+D9x+OZxJOaeRYHsGmQy7Tv9C7UxLOAFHys/dHRth77NgniHEqoTGGYMM6qjMgtz8es/x3Dw2gn8k3kOmaqLgDy/XD9ZqT0amDRFC/tW6O79DMKatOVdSkhyGGYMM6on1Go1Dl+7gN2XjiLuThxSCi/oPDQp1MYwE55wN2+KNk6tEOIbiC5eLSA3khuocqLHY5gxzKgeS8m5h1//OY7YGydwKfs8slSXAXlB+Y5qBaxk3vCyao62zq3Qs1F7tHXx4a24qNZgmDHMiDTUajWO37iI6KsncPr2WVzPu4ACWbL2s9zKqKxgLfOGj3VztHMJQA+ftmjt2pABRwbBMGOYET1SfnEx9l05g4PXTuF8xnmkFl5GsfxmuVtxAfhvwHnB06oxWjv5oZt3W3T0bMpDlKR3DDOGGVGVZRXm448rcfgz+TQS7p7H7aLLKDa6XW70JABAbQpzeMLVvBH87Jujk0cAuvsGwMbMvOYLpzqLYcYwI6oWWYV52Hflb8TeOIuEjASkFFxBkdFNnYcohTCCicoVDqbe8LVpjDbOLfCcdwD8nNx5mJKeCMOMYUakN/nFxTiUGI/YG2dwPj0BN/MvI09c0z3IBABUFrCUecDF3AfN7ZugvZsfnvNpBSdL2xqtm6SHYcYwI6pRarUa59KSsT/xNP6+HY+k7Cu4W5qMUqO08ncw+S+Zyh5KI0+4WfqghUNTdHBviWcb+sHGnNfD0X0MM4YZUa2QVZiPQ4nncfTmefyTcQG38hORq74OGGfr7C+EEeQqRyiN3eFm4YVmDo3QzqU5nvX2g4MF/z3XNwwzhhlRrZZ4Nw0x187hdGoCLmdewu3CJBTgJmBUWOE0RipbWBm5w8W8IXxtGyHAqQk6efqhsYMzb7xcRzHMGGZEkiOEwD/pNxCbHI+zty/iStZVpBUkIx8pgDyn4glVFjCDGxxMPeBh3RDN7H3Q1rUJnvFoAqWZRc2tAFU7hhnDjKjOEELg6t07iE2Ox5nbl3A16ypSC64hV30TKqN7FZ6TE0IGI7UdrIxc4GTmAS8bLzR39EEb5yZo4+YDcxM+aaC2Y5gxzIjqhcyCXBy5/g9Op1zAxbuJuJmXjHvFt1Aku/3IQ5b3z805wEp+P+g8rb3QzMEbAc6+aOvmw5sy1xIMM4YZUb2mVqtxOeM2jt+4gHN3ruBqZhJSC24gpzQFJUZpum/l9V9CyGCksoOFUQPYK1zhZumORrYN4efkg7ZujeBp04Dn6GoIw4xhRkQVKFWpEH/nOk7evIyEjCtIyrqG2wXXkVN6G8Wy9EcGHQBArYCJcIS13BmOZq7wsPJAIzsvtHTyRmtXHzhaWtXMitQDDDOGGRE9AbVajct3U3Dq1hX8k56Eq5nJSM27iXslqSgUdwDjrMfPRGUFU9jDWu4ERzNnuFm5wcfOA80dG6K1izdcrR24Z1dJDDOGGRHpwb2CPMTduopzdxJx5W4ybuTeQHrhLWSXpqFEdgcwKn78TNSmMBb2sDRyhK2pM5wtXOFp7QpfO0/4OXnBz8kDFqam+l8ZCWCYMcyIqIYJIXAr5y7OpCbinzvJSMy8gZu5t5BReBs5pXdQLMsA5LmVmI8RZCprmMIOVsYOsFM0gLOFEzxsXOFj64ZmDTzh18ADlqZ1/6bODDOGGRHVQum5OTiXlox/7lxDYuZN3My9hbSCVGSXpKFAZPz3UgMdTynQRWUJE9jBUm4PGxNHOJo5wdXKGQ1tXOFr74bmjp7wVNpDLpfuTZ4ZZgwzIpKgUlUprt67jYS067h87yauZ6UgJS8VGYV3kFOSjkJxFyqjrMcPUvkvoTaBkdoGCpkSFnI72Jo6wMHcES6WDeBh4wwfO1c0tneFt50TTOQmel67qqvK97lxDdVERESPYSw3RlNHdzR1dK+wj0qlRuK9O7iYcRNX7t7E9exUpOSmIqMoDdnF6chX30UJMgF5PmRGJRBGGShEBgoB3C0GrhYDyAJw63/zvH+BuRWMoYS5kS2sje1hq3BAA4sGcLF0hKeNC7xsXdDI3gWu1spaubfHPTMiojoouygPlzNScDnjFq5l3cat7Nu4nZ+Ou4XpyCm9iwJVJkqQCSHPrfAuKroItTGM1NYwgQ3M5LawNraFUmEHR3NHOFs4wM2mAdytneCpdEJDpSOszUyfePQmDzMyzIiIKiW/uBhX797G1bupuJaVilu5aUjLu4OMwgxkFWcgX3UPRSITKll25UZrPkAIGTyMg7F75OdPVBsPMz4ttQo4twXw7ADYeRu6GiIivbEwNYW/iyf8XTwf2zerMBdX7qYi6d5t3Mi+g1s5d5CWn457hXeRVXIP+aWZKBJZKEUOIM+DTCZgJq+ZW4MxzHQ5/S3wy5T7/7+wEhdJEhHVA0ozK7Rza4x2bo0f27d4SSNkFd1F4eCQGqgMqH1n8WqDpMOGroCISNJM89PRQKWG582jNbI8hplOvNUMEVH1qJnvU4aZLrxvGhFR9aih71OGmU4MMyKi6sEwMxzumRERVQ/umdUSVw8augIiImm5fvyBNwwzA3pg4/880XBlEBFJ0Xcv/e//uWdmQA9ufFHJO1gTEdF9RdkPvKkjYRYZGQkfHx+YmZkhMDAQMTExj+x/8OBBBAYGwszMDL6+vvjyyy/1XaIOD2z8vDQDLJ+IiKpCr2G2adMmTJ06FXPnzsXp06fRtWtX9O7dG8nJyTr7JyYmok+fPujatStOnz6NOXPmYPLkydiyZYs+yyzvwT8kVMVAQWbNLp+ISKry72q/rwuHGZcuXYpx48Zh/PjxaNGiBZYtWwZPT0+sXLlSZ/8vv/wSDRs2xLJly9CiRQuMHz8eY8eOxccff6zPMh8vJc6wyycikooLux5qqJkw09u9GYuLi3Hy5EnMmjVLqz00NBSxsbE6pzly5AhCQ0O12sLCwrBmzRqUlJTAxKT8w+OKiopQVFSkeZ+dnV2uT5XJHsr49QOA+feA4lwgbiNwfiuQceX+XpuqBECdevAAkURV8KWpc89AR9tT9auob2X7VWWe1bw+TzVPGWDnBTTqAbR8EbDzAX7+VyXmXf30Fmbp6elQqVRwdnbWand2dkZqaqrOaVJTU3X2Ly0tRXp6OlxdXctNExERgXfffbf6Cgeg84e4yA4wMgHUlXvCKxFRvZBzC0g+Auz/oIIOEg+zMg8/lE0I8cgHtenqr6u9zOzZszF9+nTN++zsbHh6Pv5RBk9EXQI4NgU6vA54dgRMLQEjefk9OSKqeTofzVjBUZOn7lsbaqgFfdWlwO1zQPwO4Or+++8fVkP3oNBbmDk6OkIul5fbC0tLSyu391XGxcVFZ39jY2M4ODjonEahUEChUFRP0WV0BedLa+4/30zpyTuEEBGVcW8HtBt1f+BH4kFg85iHOkh8AIipqSkCAwMRHR2t1R4dHY3OnTvrnCYoKKhc/z179qB9+/Y6z5fpz0Mbv9scoNVgwLYhg4yISBcL+/vnzbrNMcji9Xp8bPr06fj666+xdu1aJCQkYNq0aUhOTsaECRMA3D9EOGrUKE3/CRMm4Nq1a5g+fToSEhKwdu1arFmzBjNmzNBnmeU9HFhdp+vuR0RE2p576Pu6hk7D6PWc2bBhw5CRkYFFixYhJSUF/v7+2LVrF7y8vAAAKSkpWtec+fj4YNeuXZg2bRpWrFgBNzc3LF++HC+99FJFi9CTh8JMXpN7hUREEmYk134v9dGMZd566y289dZbOj9bt25dubbg4GCcOnVKz1VVgYmFoSsgIpIWM1ugMPP+/9fQLQE5DE+XB/+SGLnVcHUQEUlRv0//9/9qhpkBPRBmPMRIRFQ1WjdrV9XIIhlmumgd4+XoRSKiqnnge1PNMKsdmGVERFXz4AhGnjMzoAf/knBobLg6iIikqHHI//6fhxkNqGzjtx8HmCkNWwsRkdSYWgId719PzMOMhlS28W3cDFsHEZFUyU3v/5eHGQ2o7IaavIEwEdGTKbt4mntmBlR2mPHhK9mJiKhyZP/9/uQ5MwMq+0tCxjAjInoiZTsDPMxoQNwzIyJ6OmWnaXiY0YC4Z0ZE9HRq+DCj3m80LEku/kBRNqD0MHQlRETS1LAj8Ow0wK1djSxOJkRVngle+2VnZ0OpVCIrKws2NjaGLoeIiJ5QVb7PeZiRiIgkj2FGRESSxzAjIiLJY5gREZHkMcyIiEjyGGZERCR5DDMiIpI8hhkREUkew4yIiCSPYUZERJLHMCMiIsljmBERkeQxzIiISPIYZkREJHkMMyIikjyGGRERSR7DjIiIJI9hRkREkscwIyIiyWOYERGR5DHMiIhI8hhmREQkeQwzIiKSPIYZERFJHsOMiIgkT69hdu/ePYSHh0OpVEKpVCI8PByZmZkV9i8pKcE777yDVq1awdLSEm5ubhg1ahRu3bqlzzKJiEji9Bpmw4cPR1xcHHbv3o3du3cjLi4O4eHhFfbPz8/HqVOn8H//9384deoUtm7diosXL+KFF17QZ5lERCRxMiGE0MeMExIS4Ofnh6NHj6Jjx44AgKNHjyIoKAj//PMPmjVrVqn5/PXXX+jQoQOuXbuGhg0blvu8qKgIRUVFmvfZ2dnw9PREVlYWbGxsqmdliIioxmVnZ0OpVFbq+1xve2ZHjhyBUqnUBBkAdOrUCUqlErGxsZWeT1ZWFmQyGWxtbXV+HhERoTmMqVQq4enp+bSlExGRxOgtzFJTU+Hk5FSu3cnJCampqZWaR2FhIWbNmoXhw4dXmMqzZ89GVlaW5nX9+vWnqpuIiKSnymG2cOFCyGSyR75OnDgBAJDJZOWmF0LobH9YSUkJXn75ZajVakRGRlbYT6FQwMbGRutFRET1i3FVJ5g4cSJefvnlR/bx9vbGmTNncPv27XKf3blzB87Ozo+cvqSkBEOHDkViYiL27dvHgCIiokeqcpg5OjrC0dHxsf2CgoKQlZWF48ePo0OHDgCAY8eOISsrC507d65wurIgu3TpEvbv3w8HB4eqlkhERPWM3s6ZtWjRAs8//zxee+01HD16FEePHsVrr72Gfv36aY1kbN68ObZt2wYAKC0txeDBg3HixAls2LABKpUKqampSE1NRXFxsb5KJSIiidPrdWYbNmxAq1atEBoaitDQUAQEBODbb7/V6nPhwgVkZWUBAG7cuIEdO3bgxo0baNOmDVxdXTWvqoyAJCKi+kVv15kZSlWuSyAiotqrVlxnRkREVFMYZkREJHkMMyIikjyGGRERSR7DjIiIJI9hRkREkscwIyIiyWOYERGR5DHMiIhI8hhmREQkeQwzIiKSPIYZERFJHsOMiIgkj2FGRESSxzAjIiLJY5gREZHkMcyIiEjyGGZERCR5DDMiIpI8hhkREUkew4yIiCSPYUZERJLHMCMiIsljmBERkeQxzIiISPIYZkREJHkMMyIikjyGGRERSR7DjIiIJI9hRkREkscwIyIiyWOYERGR5DHMiIhI8hhmREQkeQwzIiKSPIYZERFJHsOMiIgkj2FGRESSp9cwu3fvHsLDw6FUKqFUKhEeHo7MzMxKT//GG29AJpNh2bJlequRiIikT69hNnz4cMTFxWH37t3YvXs34uLiEB4eXqlpt2/fjmPHjsHNzU2fJRIRUR1grK8ZJyQkYPfu3Th69Cg6duwIAFi9ejWCgoJw4cIFNGvWrMJpb968iYkTJ+L3339H37599VUiERHVEXrbMzty5AiUSqUmyACgU6dOUCqViI2NrXA6tVqN8PBwzJw5Ey1btnzscoqKipCdna31IiKi+kVvYZaamgonJ6dy7U5OTkhNTa1wusWLF8PY2BiTJ0+u1HIiIiI05+SUSiU8PT2fuGYiIpKmKofZwoULIZPJHvk6ceIEAEAmk5WbXgihsx0ATp48ic8++wzr1q2rsM/DZs+ejaysLM3r+vXrVV0lIiKSuCqfM5s4cSJefvnlR/bx9vbGmTNncPv27XKf3blzB87Ozjqni4mJQVpaGho2bKhpU6lU+Pe//41ly5YhKSmp3DQKhQIKhaJqK0FERHVKlcPM0dERjo6Oj+0XFBSErKwsHD9+HB06dAAAHDt2DFlZWejcubPOacLDw9GzZ0+ttrCwMISHh+PVV1+taqlERFRP6G00Y4sWLfD888/jtddew1dffQUAeP3119GvXz+tkYzNmzdHREQEXnzxRTg4OMDBwUFrPiYmJnBxcXnk6EciIqrf9Hqd2YYNG9CqVSuEhoYiNDQUAQEB+Pbbb7X6XLhwAVlZWfosg4iI6jiZEEIYuojqlJ2dDaVSiaysLNjY2Bi6HCIiekJV+T7nvRmJiEjyGGZERCR5DDMiIpI8hhkREUkew4yIiCSPYUZERJLHMCMiIsljmBERkeQxzIiISPIYZkREJHkMMyIikjyGGRERSR7DjIiIJI9hRkREkscwIyIiyWOYERGR5BkbuoDqVvas0ezsbANXQkRET6Pse7wyz5Cuc2GWk5MDAPD09DRwJUREVB1ycnKgVCof2UcmKhN5EqJWq3Hr1i1YW1tDJpM90Tyys7Ph6emJ69evP/ZR3fUBt4c2bo//4bbQxu2h7Wm3hxACOTk5cHNzg5HRo8+K1bk9MyMjI3h4eFTLvGxsbPgL+QBuD23cHv/DbaGN20Pb02yPx+2RleEAECIikjyGGRERSR7DTAeFQoEFCxZAoVAYupRagdtDG7fH/3BbaOP20FaT26PODQAhIqL6h3tmREQkeQwzIiKSPIYZERFJHsOMiIgkj2FGRESSxzDTITIyEj4+PjAzM0NgYCBiYmIMXZLeRURE4JlnnoG1tTWcnJwwcOBAXLhwQauPEAILFy6Em5sbzM3N0a1bN5w/f95AFdeciIgIyGQyTJ06VdNW37bFzZs3MXLkSDg4OMDCwgJt2rTByZMnNZ/Xp+1RWlqKefPmwcfHB+bm5vD19cWiRYugVqs1fery9jh06BD69+8PNzc3yGQybN++Xevzyqx7UVERJk2aBEdHR1haWuKFF17AjRs3nq4wQVp++OEHYWJiIlavXi3i4+PFlClThKWlpbh27ZqhS9OrsLAwERUVJc6dOyfi4uJE3759RcOGDUVubq6mz4cffiisra3Fli1bxNmzZ8WwYcOEq6uryM7ONmDl+nX8+HHh7e0tAgICxJQpUzTt9Wlb3L17V3h5eYkxY8aIY8eOicTERPHHH3+Iy5cva/rUp+3x/vvvCwcHB/Hrr7+KxMREsXnzZmFlZSWWLVum6VOXt8euXbvE3LlzxZYtWwQAsW3bNq3PK7PuEyZMEO7u7iI6OlqcOnVKdO/eXbRu3VqUlpY+cV0Ms4d06NBBTJgwQautefPmYtasWQaqyDDS0tIEAHHw4EEhhBBqtVq4uLiIDz/8UNOnsLBQKJVK8eWXXxqqTL3KyckRTZo0EdHR0SI4OFgTZvVtW7zzzjvi2WefrfDz+rY9+vbtK8aOHavVNmjQIDFy5EghRP3aHg+HWWXWPTMzU5iYmIgffvhB0+fmzZvCyMhI7N69+4lr4WHGBxQXF+PkyZMIDQ3Vag8NDUVsbKyBqjKMrKwsAIC9vT0AIDExEampqVrbRqFQIDg4uM5um3/961/o27cvevbsqdVe37bFjh070L59ewwZMgROTk5o27YtVq9erfm8vm2PZ599Fnv37sXFixcBAH///TcOHz6MPn36AKh/2+NBlVn3kydPoqSkRKuPm5sb/P39n2r71Lm75j+N9PR0qFQqODs7a7U7OzsjNTXVQFXVPCEEpk+fjmeffRb+/v4AoFl/Xdvm2rVrNV6jvv3www84deoU/vrrr3Kf1bdtcfXqVaxcuRLTp0/HnDlzcPz4cUyePBkKhQKjRo2qd9vjnXfeQVZWFpo3bw65XA6VSoUPPvgAr7zyCoD69/vxoMqse2pqKkxNTWFnZ1euz9N8zzLMdHj4OWhCiCd+NpoUTZw4EWfOnMHhw4fLfVYfts3169cxZcoU7NmzB2ZmZhX2qw/bArj/jMD27dvjP//5DwCgbdu2OH/+PFauXIlRo0Zp+tWX7bFp0yZ899132LhxI1q2bIm4uDhMnToVbm5uGD16tKZffdkeujzJuj/t9uFhxgc4OjpCLpeX++sgLS2t3F8addWkSZOwY8cO7N+/X+u5cC4uLgBQL7bNyZMnkZaWhsDAQBgbG8PY2BgHDx7E8uXLYWxsrFnf+rAtAMDV1RV+fn5abS1atEBycjKA+vW7AQAzZ87ErFmz8PLLL6NVq1YIDw/HtGnTEBERAaD+bY8HVWbdXVxcUFxcjHv37lXY50kwzB5gamqKwMBAREdHa7VHR0ejc+fOBqqqZgghMHHiRGzduhX79u2Dj4+P1uc+Pj5wcXHR2jbFxcU4ePBgnds2ISEhOHv2LOLi4jSv9u3bY8SIEYiLi4Ovr2+92RYA0KVLl3KXaVy8eBFeXl4A6tfvBgDk5+eXe+qxXC7XDM2vb9vjQZVZ98DAQJiYmGj1SUlJwblz555u+zzx0JE6qmxo/po1a0R8fLyYOnWqsLS0FElJSYYuTa/efPNNoVQqxYEDB0RKSormlZ+fr+nz4YcfCqVSKbZu3SrOnj0rXnnllToz3PhxHhzNKET92hbHjx8XxsbG4oMPPhCXLl0SGzZsEBYWFuK7777T9KlP22P06NHC3d1dMzR/69atwtHRUbz99tuaPnV5e+Tk5IjTp0+L06dPCwBi6dKl4vTp05rLlyqz7hMmTBAeHh7ijz/+EKdOnRI9evTg0Hx9WLFihfDy8hKmpqaiXbt2muHpdRkAna+oqChNH7VaLRYsWCBcXFyEQqEQzz33nDh79qzhiq5BD4dZfdsWv/zyi/D39xcKhUI0b95crFq1Suvz+rQ9srOzxZQpU0TDhg2FmZmZ8PX1FXPnzhVFRUWaPnV5e+zfv1/nd8Xo0aOFEJVb94KCAjFx4kRhb28vzM3NRb9+/URycvJT1cXnmRERkeTxnBkREUkew4yIiCSPYUZERJLHMCMiIsljmBERkeQxzIiISPIYZkREJHkMMyIikjyGGRERSR7DjIiIJI9hRkREkvf/4b5/Dk9MxpAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 480x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "eqn = eq_continuity & eq_momentum & eq_eth\n",
    "     \n",
    "viewer = fipy.viewers.matplotlibViewer.matplotlib1DViewer.Matplotlib1DViewer(vars=(dens, vel*unit_velocity/1.e6, Eth*(gamma-1)),\n",
    "                    datamin=-0.5, datamax=1.,legend='upper left', title = '1D hydro solution')\n",
    "viewer.axes.legend([r'Density ($10^{-12}\\ \\mathrm{kg\\ m^{-3}}$)',r'Velocity ($\\mathrm{10\\ km/s}$)',r'Pressure (0.03 Pa)'])\n",
    "viewer.plot()\n",
    "\n",
    "steps=50 # 5000\n",
    "\n",
    "# 60% of maximal sound speed in the domain\n",
    "cflcond = 0.6 * dx/np.max(np.sqrt(gamma**2*Eth.value/dens.value))\n",
    "timeStepDuration = cflcond\n",
    "t=np.arange(steps)*timeStepDuration\n",
    "    \n",
    "for step in np.arange(steps):\n",
    "    Eth.updateOld()\n",
    "    dens.updateOld()\n",
    "    vel.updateOld()\n",
    "\n",
    "    for sweep in range(5):\n",
    "        EthPrevious[:] = Eth\n",
    "        densPrevious[:] = dens\n",
    "        velPrevious[:] = vel\n",
    "\n",
    "        eqn.cacheMatrix()\n",
    "        residual = eqn.sweep(dt=timeStepDuration)\n",
    "\n",
    "        matrixDiagonal[:] = eqn.matrix.takeDiagonal()[mesh.numberOfCells:2 * mesh.numberOfCells]\n",
    "        Eth[:] = relaxation * Eth + (1 - relaxation) * EthPrevious\n",
    "        dens[:] = relaxation * dens + (1 - relaxation) * densPrevious\n",
    "        vel[:] = relaxation * vel + (1 - relaxation) * velPrevious\n",
    "\n",
    "    if step%100 ==0:\n",
    "        print('done step',step, ', time', step*timeStepDuration*unit_time)\n",
    "    viewer.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6c8a168-60d0-4fe8-b120-8053b375da26",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}