{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical Solution of Ordinary Differential Equations: Part II\n",
    "\n",
    "This notebook continues from part I, and discusses boundary value ODEs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import some useful stuff\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "plt.rcParams[\"figure.figsize\"] = (7,5)   # Make default figure size bigger"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Boundary Value Problems\n",
    "\n",
    "This second notebook deals with more the class of ODEs known as boundary value problems. To remind you, and ODE is an initial value problem if all of the boundary conditions are specified at the same value of the independent variable. An ODE where this is not the case is generically referred to as a boundary value problem.\n",
    "\n",
    "As a simple example, consider the ODE describing the motion of a point mass (for example a ball) traveling in a constant gravitational field. The height of the ball at any time is $z$, and gravity exerts a constant acceleration $-g$ in the $z$ direction. The position of the ball then obeys the second order ODE\n",
    "\\begin{equation}\n",
    "\\frac{d^2 z}{dt^2} = -g.\n",
    "\\end{equation}\n",
    "Since this is a second order ODE, it requires two boundary conditions. For example, one might specify that at time $t=0$ the ball is initially at zero height ($z=0$ at $t=0$) and is thrown upward at a constant speed $v$ ($dz/dt = v$ at $t=0$). Finding the subsequent height versus time is then an initial value problem.\n",
    "\n",
    "However, one could also specify the system a different way. One could say that the ball is at height $z=0$ at $t=0$, and returns to that height after $t=5$ seconds, so that $z=0$ at $t=5$. This also provides the two boundary conditions required to uniquely specify the height of the ball at all times, but mathematically and computationally the structure of the problem is now quite different. This is an exaple of a class of ODE known as a boundary value problem; it is an ODE where the solution is to be found between two different values of the independent variable, with boundary conditions specified at both of them. Nor is this the only way the boundary value problem could be specified. For example, we could say that the ball is at $z=0$ at $t=0$, and that its velocity is $dz/dt = 2$ m/s at $t=2$ s. This is also a valid set of boundary conditions, and also constitutes a boundary value problem. \n",
    "\n",
    "One could even imagine more complex cases where one of the boundary conditions isn't expressed at a fixed time. For example, we could ask that the ball reach zero velocity at a height $z = 2$ m, so that our two boundary conditions are $z = 0$ at $t = 0$ and $dz/dt = 0$ when $z = 2$ m; in this case the second boundary condition doesn't explicitly specify time at all, and instead specifies a relationship between the dependent variable and its derivative at a particular point, in the form of an algebraic equation:\n",
    "\\begin{equation}\n",
    "\\left(\\frac{dz}{dt}\\right)_{z=2\\,\\mathrm{m}} = 0.\n",
    "\\end{equation}\n",
    "Or we could specify that the ball must travel a total distance of 2 m (in either direction, up or down) in a time of one second, in which case the boundary condition is being expressed as a requirement on the integral of the solution:\n",
    "\\begin{equation}\n",
    "\\int_0^{2\\,\\textrm{s}} \\left|\\frac{dz}{dt}\\right| \\, dt = 2\\,\\textrm{m}.\n",
    "\\end{equation}\n",
    "Nonetheless, there is clearly a valid solution to the ODE that satisfies these constraints. In some cases there may be only a single valid solution, while in others there may be more than one. Although these two examples are not techincally cases where conditions are specified at the boundary, we also usually refer to problems like this as boundary value problems, since the methods use to tackle them are the same as for standard boundary value problems.\n",
    "\n",
    "There is a vast universe of methods for handling ODEs of this type, and in general they are much harder than initial value problems. A good physical understanding of the nature of the problem is often required to select a good numerical method. In this notebook we will introduce two of the most common approaches."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shooting\n",
    "\n",
    "The first method we will consider is called shooting. To motivate the idea of shooting, let us consider a case that illustrates the physical picture behind the method. An archer is firing an arrow at a target that is at a distance $x$ and at the same elevation as the archer's position. When fully drawn, the bow launches the arrow with a fixed velocity $v$. (We'll pretend it's a crossbow rather than a regular bow, so the launch velocity is determined by a cranking system and is close to the same for every shot.) The bow can be rotated freely so that the arrow is launched at an angle $\\theta$ relative to horizontal. At what angle $\\theta$ should the shot be fired so as to hit the target? To make the problem non-trivial, let us suppose that the arrow is subject not only to gravity, but to aerodynamic drag, which exerts a velocity-dependent force in the direction opposite the arrow's motion.\n",
    "\n",
    "One could imagine solving this problem in real life by firing the bow at a bunch of different angles. Suppose that $z>0$, so the target is above the location of the bow. For angles $\\theta$ close to zero or close to $90^\\circ$, the shot will travel negligible horizontal distance before falling below its initial height (and hitting the ground if the bow is close to the ground), and will therefore miss. As the angle is tipped away from purely horizontal or purely vertical the shot travels further. One could gradually change $\\theta$, noting whether the arrows overshoots or undershoots the target, adjusting the angle until the target is hit. This is exactly the idea behind shooting.\n",
    "\n",
    "Now let us formalize this idea. The ODEs describing the trajectory of the arrow are\n",
    "\\begin{eqnarray}\n",
    "\\frac{dx}{dt} & = & v_x \\\\\n",
    "\\frac{dz}{dt} & = & v_z \\\\\n",
    "\\frac{dv_x}{dt} & = & - \\frac{C_d A \\rho_{\\mathrm{air}}}{2 m} v_x \\sqrt{v_x^2+v_z^2} \\\\\n",
    "\\frac{dv_z}{dt} & = & - \\frac{C_d A \\rho_{\\mathrm{air}}}{2 m} v_z \\sqrt{v_x^2 + v_z^2} - g,\n",
    "\\end{eqnarray}\n",
    "where $m$ is the arrow mass, $g$ is the gravitational acceleration, $C_d$ is the coefficient of drag, $A$ is the cross-sectional area of the arrow, $\\rho_{\\mathrm{air}}$ is the density of air, and we have assumed a drag law such that the magnitude of the drag force scales with total velocity as $v^2$. This is a fourth-order system for the variables $x$, $z$, $v_x$, and $v_z$.\n",
    "\n",
    "For a shot with a specified angle $\\theta$, the problem is simply an initial value problem, with the four boundary conditions\n",
    "\\begin{eqnarray}\n",
    "\\left. x \\right|_{t=0} & = & 0 \\\\\n",
    "\\left. z \\right|_{t=0} & = & 0 \\\\\n",
    "\\left. v_x \\right|_{t=0} & = & v \\cos\\theta \\\\\n",
    "\\left. v_z \\right|_{t=0} & = & v \\sin\\theta \\\\\n",
    "\\end{eqnarray}\n",
    "These four boundary conditions fully specify the problem, and thus for any choice of $\\theta$ we can immediately generate the trajectory of the arrow. In paticular, we can compute the horizontal distance $x$ travelled by the arrow for any given value of $\\theta$.\n",
    "\n",
    "We can therefore break the problem up into two parts:\n",
    "\n",
    "1. Write a function $x(\\theta)$ computes the horizontal distance $x$ as a function of $\\theta$. This is an algebraic function of one variable.\n",
    "2. Solve the algebraic system $x(\\theta) = x_{\\mathrm{target}}$ to find the value of $\\theta$ that hits the target.\n",
    "\n",
    "This is the essence of the shooting method: we guess one or more parameters so that we can convert our ODE to an initial value problem. We then search for a solution as a function of those parameters, exactly the same way we would solve any non-linear system of equations (e.g., using Newton's method, the secant method, or the like).\n",
    "\n",
    "One thing to note before proceeding: if the system of ODEs we are solving is linear (meaning that the dependent variable appears only linearly, i.e., there are terms involving $x$ and $v_x$, but none involving $e^x$ or $v_x^2$ or anything like that), then a significant simplification is possible. Since this occurs rarely in practice, and is not the case even for the simple example we are examining, we will not bother with this reduced case. If you encounter a problem of this sort, known as a linear boundary value problem, you can find standard solutions for it fairly easily."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Writing the Objective Function\n",
    "\n",
    "Let us begin with the first step: writing a function that, for a given choice of $\\theta$, returns the horizontal distance travelled by the arrow before it hits the ground. More generally, we can think of this step as writing the \"objective function\": the function that represents the condition we're trying to satisfy, in this case that the horizontal distance travelled by the arrow matches the distance to the target. This in itself requires some thought, because we have to watch the value of the solution, and terminate when the height becomes negative. We will then need to go back and refine the solution to identify the exact horizontal distance at which the arrow goes from $z>0$ to $z<0$. We will accomplish this by means of a loop where we advance the solution in time until the height is negative, then go back and re-integrate with a smaller time step, until the time step becomes small enough or the height close enough to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set the drag coefficient term and g; note that the value used is a realistic one for an arrow; units are SI\n",
    "drag_fac = 0.003\n",
    "g = 9.8\n",
    "\n",
    "# First define the derivatives function\n",
    "def f(t, q):\n",
    "    # t = time, q[0] = x, q[1] = z, q[2] = vx, q[3] = vz\n",
    "    x = q[0]\n",
    "    z = q[1]\n",
    "    vx = q[2]\n",
    "    vz = q[3]\n",
    "    # Derivatives of position\n",
    "    dqdt = np.zeros(4)\n",
    "    dqdt[0] = vx\n",
    "    dqdt[1] = vz\n",
    "    # Derivatives of velocity\n",
    "    dqdt[2] = -drag_fac * vx * np.sqrt(vx**2+vz**2)\n",
    "    dqdt[3] = -drag_fac * vz * np.sqrt(vx**2+vz**2) - g\n",
    "    return dqdt\n",
    "\n",
    "# Set up an integrator using this function; this time use Dormand-Price 4th/5th order\n",
    "from scipy.integrate import ode\n",
    "integrator = ode(f).set_integrator('dopri5')\n",
    "\n",
    "# Now define a function that takes theta as an argument, and returns the horizontal distance;\n",
    "# the optional arguments are as follows:\n",
    "# v = launch velocity of arrow; a reasonable value for a crossbow\n",
    "# dt = guess for initial time step\n",
    "# dt_lim = minimum time step\n",
    "# z_lim = accuracy with which we want the height to go to zero\n",
    "# verbose = print out z as the calculation runs\n",
    "def hdist(theta, v=100.0, dt=10.0, dt_lim=1.0e-10, z_lim=1.0e-6, verbose=False):\n",
    "    \n",
    "    # Set initial conditions\n",
    "    q = np.array([0.0, 0.0, v*np.cos(theta), v*np.sin(theta)])\n",
    "    t0 = 0.0\n",
    "\n",
    "    # Set the initial conditions and parameters\n",
    "    integrator.set_initial_value(q, t0)\n",
    "    \n",
    "    # Loop\n",
    "    while True:\n",
    "        \n",
    "        # Save previous time step\n",
    "        tlast = integrator.t\n",
    "        qlast = q\n",
    "\n",
    "        # Take a step\n",
    "        q = integrator.integrate(integrator.t+dt)\n",
    "        if verbose:\n",
    "            print(\n",
    "                \"z = {:e}, dt = {:e}, t = {:f}, x = {:f}, vx = {:f}, vz = {:f}\".\n",
    "                format(q[1], dt, integrator.t, q[0], q[2], q[3]))\n",
    "\n",
    "        # Check termination condition on height and velocity\n",
    "        if abs(q[1]) < z_lim and q[3] < 0.0:\n",
    "            break\n",
    "            \n",
    "        # If we overshot and the vertical position is negative, rewind to previous\n",
    "        # state, and reduce time step\n",
    "        if q[1] < 0.0:\n",
    "            dt = dt/2.0\n",
    "            integrator.set_initial_value(qlast, tlast)\n",
    "            q = qlast\n",
    "\n",
    "        # Check termination condition on dt\n",
    "        if dt < dt_lim:\n",
    "            break\n",
    "            \n",
    "    # Return horizontal distance\n",
    "    return q[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "z = 6.663472e+00, dt = 1.000000e+01, t = 10.000000, x = 349.075623, vx = 18.429464, vz = -39.864262\n",
      "z = -5.076748e+02, dt = 1.000000e+01, t = 20.000000, x = 444.593351, vx = 3.807839, vz = -56.280967\n",
      "z = -2.324265e+02, dt = 5.000000e+00, t = 15.000000, x = 414.814755, vx = 8.740223, vz = -53.030754\n",
      "z = -1.048640e+02, dt = 2.500000e+00, t = 12.500000, x = 388.014684, vx = 12.924565, vz = -48.501445\n",
      "z = -4.641440e+01, dt = 1.250000e+00, t = 11.250000, x = 370.268319, vx = 15.524557, vz = -44.825347\n",
      "z = -1.910104e+01, dt = 6.250000e-01, t = 10.625000, x = 360.126115, vx = 16.942708, vz = -42.521077\n",
      "z = -6.011266e+00, dt = 3.125000e-01, t = 10.312500, x = 354.717040, vx = 17.678138, vz = -41.238509\n",
      "z = 3.797822e-01, dt = 1.562500e-01, t = 10.156250, x = 351.925681, vx = 18.051901, vz = -40.563049\n",
      "z = -6.011266e+00, dt = 1.562500e-01, t = 10.312500, x = 354.717040, vx = 17.678138, vz = -41.238509\n",
      "z = -2.802549e+00, dt = 7.812500e-02, t = 10.234375, x = 353.328661, vx = 17.864534, vz = -40.903671\n",
      "z = -1.208057e+00, dt = 3.906250e-02, t = 10.195312, x = 352.629001, vx = 17.958097, vz = -40.734086\n",
      "z = -4.133024e-01, dt = 1.953125e-02, t = 10.175781, x = 352.277799, vx = 18.004970, vz = -40.648750\n",
      "z = -1.655088e-02, dt = 9.765625e-03, t = 10.166016, x = 352.101855, vx = 18.028428, vz = -40.605945\n",
      "z = 1.816680e-01, dt = 4.882812e-03, t = 10.161133, x = 352.013796, vx = 18.040163, vz = -40.584509\n",
      "z = -1.655088e-02, dt = 4.882812e-03, t = 10.166016, x = 352.101855, vx = 18.028428, vz = -40.605945\n",
      "z = 8.257165e-02, dt = 2.441406e-03, t = 10.163574, x = 352.057833, vx = 18.034295, vz = -40.595230\n",
      "z = -1.655088e-02, dt = 2.441406e-03, t = 10.166016, x = 352.101855, vx = 18.028428, vz = -40.605945\n",
      "z = 3.301366e-02, dt = 1.220703e-03, t = 10.164795, x = 352.079845, vx = 18.031361, vz = -40.600588\n",
      "z = -1.655088e-02, dt = 1.220703e-03, t = 10.166016, x = 352.101855, vx = 18.028428, vz = -40.605945\n",
      "z = 8.232206e-03, dt = 6.103516e-04, t = 10.165405, x = 352.090850, vx = 18.029895, vz = -40.603267\n",
      "z = -1.655088e-02, dt = 6.103516e-04, t = 10.166016, x = 352.101855, vx = 18.028428, vz = -40.605945\n",
      "z = -4.159132e-03, dt = 3.051758e-04, t = 10.165710, x = 352.096353, vx = 18.029161, vz = -40.604606\n",
      "z = 2.036589e-03, dt = 1.525879e-04, t = 10.165558, x = 352.093602, vx = 18.029528, vz = -40.603936\n",
      "z = -4.159132e-03, dt = 1.525879e-04, t = 10.165710, x = 352.096353, vx = 18.029161, vz = -40.604606\n",
      "z = -1.061259e-03, dt = 7.629395e-05, t = 10.165634, x = 352.094977, vx = 18.029345, vz = -40.604271\n",
      "z = 4.876681e-04, dt = 3.814697e-05, t = 10.165596, x = 352.094289, vx = 18.029436, vz = -40.604104\n",
      "z = -1.061259e-03, dt = 3.814697e-05, t = 10.165634, x = 352.094977, vx = 18.029345, vz = -40.604271\n",
      "z = -2.867945e-04, dt = 1.907349e-05, t = 10.165615, x = 352.094633, vx = 18.029391, vz = -40.604187\n",
      "z = 1.004370e-04, dt = 9.536743e-06, t = 10.165606, x = 352.094461, vx = 18.029413, vz = -40.604146\n",
      "z = -2.867945e-04, dt = 9.536743e-06, t = 10.165615, x = 352.094633, vx = 18.029391, vz = -40.604187\n",
      "z = -9.317871e-05, dt = 4.768372e-06, t = 10.165610, x = 352.094547, vx = 18.029402, vz = -40.604167\n",
      "z = 3.629152e-06, dt = 2.384186e-06, t = 10.165608, x = 352.094504, vx = 18.029408, vz = -40.604156\n",
      "z = -9.317871e-05, dt = 2.384186e-06, t = 10.165610, x = 352.094547, vx = 18.029402, vz = -40.604167\n",
      "z = -4.477478e-05, dt = 1.192093e-06, t = 10.165609, x = 352.094526, vx = 18.029405, vz = -40.604161\n",
      "z = -2.057281e-05, dt = 5.960464e-07, t = 10.165609, x = 352.094515, vx = 18.029406, vz = -40.604159\n",
      "z = -8.471829e-06, dt = 2.980232e-07, t = 10.165608, x = 352.094510, vx = 18.029407, vz = -40.604157\n",
      "z = -2.421339e-06, dt = 1.490116e-07, t = 10.165608, x = 352.094507, vx = 18.029407, vz = -40.604157\n",
      "z = 6.039069e-07, dt = 7.450581e-08, t = 10.165608, x = 352.094506, vx = 18.029408, vz = -40.604156\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "352.09450559489915"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Demonstrate the method using the default arguments\n",
    "hdist(np.pi/4, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x10d4e9a50>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcAAAAFECAYAAACwBO69AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VfX9x/HXJ5tASBhhJWzCBhEjgjhx79HW0dq6Wtu6\ntbba1rZ2j19bbeto3Wite+EWERcqe8gm7BnCSAgkZH5+f+RiU8TkArk59+a+n4/HeST35Jx737mP\nC598v+d7vl9zd0REROJNQtABREREgqACKCIicUkFUERE4pIKoIiIxCUVQBERiUsqgCIiEpdUAEVE\nJC6pAIqISFxSARQRkbiUFHSAg9GxY0fv1atX0DFERCSKzJw5c4u7Zzd2XEwXwF69ejFjxoygY4iI\nSBQxs9XhHKcuUBERiUsqgCIiEpdUAEVEJC6pAIqISFxSARQRkbikAigiInFJBVBEROKSCqCIiMQl\nFUAREYlLMT0TjEhzK6+sYeuuCkrKqygpq2LH7ioqa5zqmlqqamoBSElKICUxkZSkBNqmJdG+dQpZ\n6Sm0S08mKVF/c4pECxVAkb3srqphedFOlhaWsrRwJyuKdrKheDfri8vZtqvygJ83waBz2zRyslqR\n064VvTu2Jq9TBv07t6FXx9YkqziKNCsVQIl764vLmbZyK7PXFDN7TTGLNu6gutYBSEowenZIJ7dd\nOkNzMsnJSiM7I5XMVilkpSeTkZZEalIiyYn2eeuusrqWyupaKqprKCmvorisiuKySopKK1hfvJv1\nxWXMXL2dCXM34HUvQ0piAoO6tWVEbiYjemSR37M93dunB/WWiMQFFUCJO7uravho2RY+XFbEhwVb\nWFG0C4DWKYkc0j2Lq47pw+BubRnQOSOiLbPyyrqW5rLNpSzaWMrctcU8O3Md4z+pm8c3J6sVR/bt\nwJH9OnBMXjYd2qRGJIdIvDLf8ydoDMrPz3etBiHh2FVRzbuLN/Pm/E1MXrKZssoaWiUnckSf9hyd\nl82YPh0Y0CWDxAQLNGdNrbO0sJTpq7bxccFWPlmxlZLyKhIMRvZoxwmDOnPykM70zW4TaE6RaGZm\nM909v9HjVAClpaqtdaau3MazM9fyxmebKK+qoWObVE4Z0plTh3ZhVO/2pCYlBh2zQbW1zvwNJUxa\ntJl3FhWyYMMOAAZ3bcvZI7px5vCu5LZTV6lIfYEXQDNLAz4AUqnran3O3X9hZo8CxwIloUMvc/c5\nZmbA34DTgbLQ/lkNvYYKoOxLcVklT05by3+mrWbttnIyUpM485BunHdoDof1bBd4K+9gbCgu5835\nm3hl3gZmrykGYHSf9lx0eA9OHdqFtOToLugizSEaCqABrd19p5klAx8BNwDfA1519+f2Ov504Drq\nCuARwN/c/YiGXkMFUOpbVljKw1NW8eLsdeyuqv28MJwypAutUlpeYVi7rYyX56zn2ZnrWL21jLZp\nSZx3aA7fOrKXukglroVbACM2CMbrKuvO0MPk0NZQtT0HeCx03qdmlmVmXd19Y6QySsswf30J/3h3\nGW8tKCQ1KYHzDs3hsrG9GNilbdDRIqp7+3SuHZfH1cf149OVW3l6+lqenLaW8Z+s5rgB2Vx2ZC+O\nycsmIYZbvCKRFNFrgGaWCMwE+gH3uPutoS7QMUAFMAm4zd0rzOxV4A/u/lHo3EnAre4+Y6/nvAq4\nCqBHjx6HrV69OmL5JbrNX1/CnROXMmnxZjLSkrh8bG8uO7IX7VunBB0tMEWlFfxn6hr+PXU1RaUV\n9O/chquP68eZw7vqJnyJG4F3ge4VJgt4kbouzq3AJiAFuB9Y7u6/MrPXgN/vVQB/5O4zv+x51QUa\nn9ZuK+PPby/h5TkbyEpP5sqxvbl0bC/apiUHHS1qVFbX8uq8Dfzz/eUsLdxJj/bpfO/Yvnz1sFxS\nklQIpWULvAu0PncvNrP3gFPd/c+h3RVm9ghwS+jxOqB7vdNygQ3NkU9iw47dVfz9nWU89slqzOCa\n4/vy3WP7qvDtQ0pSAuePzOXcETm8s6iQeyYX8JMXP+Pe9wq44YQ8zjs0Ry1CiXsRK4Bmlg1UhYpf\nK+BE4I97ruuFBsmcC8wPnTIBuNbMnqJuEEyJrv8JgLvz8pwN/Oa1RWzdVcHXDsvlppP60zWzVdDR\nol5CgnHykC6cNLgz7y8t4i9vL+WHz83jn+8v5wcnD+C0oV2o+6coEn8i2QLsCowPXQdMAJ5x91fN\n7N1QcTRgDnWjQgFep24EaAF1t0FcHsFsEiMKNpdy+0vz+XTFNg7pnsUjlx3OsNzMoGPFHDPjuAGd\nOLZ/Nm8t2MRf3l7K1U/MIr9nO24/czAjumcFHVGk2elGeIlK1TW1/OuDFfztnWW0Sknk1lMHctHh\n3TWisYnU1DrPzFjLX95ewpadlZw7ohu3njZQrWppEaJqEEykqAC2TEsLS/nhs3OZu66E04Z24dfn\nDqWj5sGMiJ0V1dz3XgEPfLiSpATjxhPzuHxsb61MITFNBVBijrvz0Ecr+dObS2iTlsSvzhnCmcO7\nBR0rLqzZWsYvX1nApMWb6d+5Db8+ZyhH9OkQdCyRA6ICKDFly84KfvjsXCYvKeLEQZ35w1eGqdUX\ngIkLC7ljwgLWF5dz8age/Pj0gRplKzEnqm6DEGnIxwVbuOHpOZSUV/Grc4bwzdE9NTIxICcN7sxR\n/Tpy1ztLeeDDFby7uJDfnDuMkwZ3DjqaSJNTR78Ext25970CLnloKm3Tknjp6rF8a0wvFb+AtUpJ\n5MenD+Kla8bSLj2F7zw2g+ufnE1xWWXQ0USalAqgBGJnRTXX/GcWf3pzCacN68qEa49icLeWPXdn\nrBmem8WEa4/iphP78/pnGznlrg/4YGlR0LFEmowKoDS7VVt2cd49U3hz/iZ+cvpA7r74UFqnqjc+\nGqUkJXDDiXm8ePVYMtKS+dbD0/jZS/Mpr6wJOprIQVMBlGY1beU2zr13Clt2VvD4lUdw1TF91eUZ\nA4blZvLqdUdxxdjePP7pas655yOWFpYGHUvkoKgASrN5ec56LnlwKu1bp/DSNWMZ269j0JFkP6Ql\nJ/Lzswbz2BWj2LarkrP+8RFPTF1NLI8kl/imAigR5+78Y9IybnhqDof2yOKF7x9Jzw6tg44lB+iY\n/tm8ccMxjOrdnp++OJ9r/zObnRXVQccS2W8qgBJRtbXOLyYs4C8Tl3LeoTk8duUostLjd72+liI7\nI5Xxl4/i1lMH8uaCTZx9t7pEJfaoAErEVNXUctMzc3jsk9VcdUwf/nrBIaQmJQYdS5pIQoLx/eP6\n8sS3j2BHeTXn3D2Fl+esDzqWSNhUACUidlfV8N3HZ/LynA388JQB/Pi0gRrs0kKN7tOB168/imE5\nmdzw1Bx+9cpCqmtqg44l0igVQGlyZZXVXP7IdCYv2cxvzh3KNcf3U/Fr4Tq1TeOJ7xzB5WN78fCU\nlVz2yHTdOC9RTwVQmlRZZTVXPjqDqSu3cucFI7hkdM+gI0kzSU5M4BdnDeFPXx3OtJXbOPvuKbou\nKFFNBVCaTHllzefF768XjODcQ3OCjiQBuCC/O099dzTlVTWcf+/HvLdkc9CRRPZJBVCaRHllDVc8\nOl3FTwAY2aMdr1x7FD07pHPFo9N5/JNVQUcS+QIVQDloldW1fP+JmXy6cit/ueAQFT8BoEtmGs98\ndwzjBnbiZy8v4FevLKSmVjfNS/RQAZSDUlPr/ODZuby3pIjfnTeM8w7NDTqSRJHWqUn865v5XDG2\nNw9PWcn3/z2T3VWaR1SigwqgHDB35xcT5vPK3A3ceupALh7VI+hIEoUSE4yfnzWYO84azMRFhVzy\n4FSNEJWooAIoB+yvE5fy70/X8N1j+/D94/oGHUei3GVje3P3xSOZt66Er/7zE9YXlwcdSeKcCqAc\nkCenreEf7xZwYX53bjt1YNBxJEacMbwr468YRWHJbs6/dwrLdJuEBEgFUPbb+0uLuP2l+RzbP5vf\nnjdUN7nLfhnTtwPPfG8MtQ4X/OsT5q0rDjqSxCkVQNkvizbu4JonZtG/cwb3fGMkSYn6CMn+G9S1\nLc99bwxt0pL4+gNT+WT51qAjSRzS/14StsIdu7ni0em0SU3i4cvyaaNV3OUg9OzQmme/eyRdM9O4\n9JFpTFpUGHQkiTMRK4BmlmZm08xsrpktMLNfhvb3NrOpZrbMzJ42s5TQ/tTQ44LQz3tFKpvsv91V\nNVz12Ax2lFfx8GWH0zWzVdCRpAXokpnG098dw8AuGXzv3zN5c/6moCNJHIlkC7ACGOfuhwAjgFPN\nbDTwR+BOd88DtgNXho6/Etju7v2AO0PHSRRwd3764nzmrivhzgtHMLhb26AjSQvSvnUK//72EQzL\nyeSa/8zi1Xkbgo4kcSJiBdDr7Aw9TA5tDowDngvtHw+cG/r+nNBjQj8/wTS6Iio8MmUVz89ax40n\n5nHykC5Bx5EWqG1aMo9deQQje2Rx/ZOzeWm21hWUyIvoNUAzSzSzOcBmYCKwHCh29+rQIeuAPfNm\n5QBrAUI/LwE67OM5rzKzGWY2o6ioKJLxBZhSsIXfvr6Ikwd35vpxeUHHkRasTWoS468YxRG9O3DT\nM3N4Yda6oCNJCxfRAujuNe4+AsgFRgGD9nVY6Ou+WntfmDjQ3e9393x3z8/Ozm66sPIF64vLufY/\ns+jTsTV/vXAECQlqkEtkpack8fBlhzOmTwdueXauVpiXiGqWUaDuXgy8B4wGssxsz/DBXGBPh/86\noDtA6OeZwLbmyCdfVFldy7X/mUVVjXP/tzTiU5pPq5REHrw0n8N7teemp+fomqBETCRHgWabWVbo\n+1bAicAiYDLw1dBhlwIvh76fEHpM6Ofvurumjg/IH99czOw1xfzxK8Pp3bF10HEkzuxpCR7Wsx03\nPDWHN+dvDDqStECRbAF2BSab2TxgOjDR3V8FbgVuNrMC6q7xPRQ6/iGgQ2j/zcBtEcwmDXhz/iYe\n+mgll47pyRnDuwYdR+JU69QkHrl8FIfkZnLdk7OZrIV1pYlZLDey8vPzfcaMGUHHaFHWbC3jjH98\nSO+OrXn2e2NITUoMOpLEuZLyKr7+wKcUbN7J+CtGMbrPF8bGifwPM5vp7vmNHaeZYORzVTW1XPfU\nbAy45+sjVfwkKmS2SuaxK0aR264V3x4/g7lrNXeoNA0VQPnc3yctY+7aYn5//nC6t08POo7I5zq0\nSeWJb4+mXetkLn1kGku1ioQ0ARVAAWD6qm3cM7mAr4zM1XU/iUpdMtN44srRpCQm8K2HprFue1nQ\nkSTGqQAKpburuOnpOeS0a8UdZw8OOo7Il+rRIZ3HrhxFWWU133p4Gtt2aWV5OXAqgMIvJixgQ3E5\nd104goy05KDjiDRoYJe2PHTZ4azfXs7lj0xjV0V14yeJ7IMKYJx7c/5GXpi1nmvH5XFYz/ZBxxEJ\ny+G92nPP10cyf8MOvvfvmVRW1wYdSWKQCmAc27arkttfms/QnLZcN65f0HFE9suJgzvz+/OH8eGy\nLdz2wjxi+ZYuCYbmt4pjd0xYQEl5FY9feQTJWtldYtAF+d3ZWLybO99ZSm5WK24+eUDQkSSGqADG\nqTfnb2LC3A3cdGJ/BnXV+n4Su64/oR8bisv5+7sFdMtqxUWjegQdSWKECmAc2h7q+hzctS1XH983\n6DgiB8XM+M15Q9m4Yzc/fWk+nTPTOH5Ap6BjSQxQv1cc+tWrCykuq+T/vjZcXZ/SIiQnJnDvN0Yy\noHMG1z4xi0UbdwQdSWKA/veLMx8uK+LF2ev5/nF9GdItM+g4Ik2mTWoSD12WT5u0JK58dDqbS3cH\nHUminApgHNldVcPPXppPrw7pXHO8Rn1Ky9M1sxUPXXo428uq+M5jMymvrAk6kkQxFcA4cu/kAlZt\nLeM35w4jLVkTXUvLNDQnk7suGsG8dcX84Nk51Nbq9gjZNxXAOFGweSf3vb+cc0d046i8jkHHEYmo\nU4Z04cenDeT1zzZx1ztLg44jUUqjQOOAu/PTFz+jVXIiPz1Dc31KfPjO0X1YVriTv79bQP8uGZw5\nvFvQkSTKqAUYB16cvZ6pK7dx22mDyM5IDTqOSLPYc3vEYT3bccuzc5m/viToSBJlVABbuJ0V1fz+\njcUckpvJRYd3DzqOSLNKTUrkn5ccRrv0FL7z2AyKSiuCjiRRRAWwhbv73QKKSiu44+whJCRY0HFE\nml12RioPfCuf7WWVmjhb/keD1wDN7PwwnmO3u7/eRHmkCa3asouHP1rJ+SNzOLRHu6DjiARmaE4m\n//fVQ7juydn86tUF/ObcYUFHkijQ2CCYB4CXgYaaDscAKoBR6DevLSQ50bjt1IFBRxEJ3FmHdGP+\n+hL+9cEKhuVkcuHhmjM03jVWAN9w9ysaOsDM/t2EeaSJvL+0iHcWbebWUwfSqW1a0HFEosIPTxnA\ngg07+NlLCxjQpS0jumcFHUkC1OA1QHe/pLEnCOcYaV5VNbX86pUF9OqQzhVH9Qo6jkjUSEpM4B8X\nH0qntql87/GZGhQT58IaBGNmiWZ2tpldb2Y379kiHU4OzNPT17K8aBc/OX0QqUma8UWkvnatU/jX\nNw+juLyS65+cTXWNBsXEq3BHgb4CXAZ0ADLqbRJldlVUc9c7yzi8VztOGtw56DgiUWlIt0x+c+4w\nPlmxlb9M1Ewx8SrcmWBy3X34/jyxmXUHHgO6ALXA/e7+NzO7A/gOUBQ69Cd7RpGa2Y+BK4Ea4Hp3\nf2t/XlPggQ9XsGVnBfd/6zDMdNuDyJf56mG5zFy9nfveW87IHvqDMR6F2wJ8w8xO3s/nrgZ+4O6D\ngNHANWa2Zx6uO919RGjbU/wGAxcBQ4BTgXvNTP13+6GotIL7P1jB6cO6MFK3PYg06hdnDWZoTltu\nfmYOq7fuCjqONLNwC+CnwItmVm5mO8ys1MwaXHHS3Te6+6zQ96XAIiCngVPOAZ5y9wp3XwkUAKPC\nzCfA3yYtpbK6lh+eotseRMKRlpzIfd84jAQzvvfvWeyu0vJJ8STcAvgXYAyQ7u5t3T3D3duG+yJm\n1gs4FJga2nWtmc0zs4fNbE9TJQdYW++0deyjYJrZVWY2w8xmFBUV7f3juLW8aCdPTlvL14/oQe+O\nrYOOIxIzurdP584LD2HRxh386tWFQceRZhRuAVwGzHf3/V5Yy8zaAM8DN7r7DuA+oC8wAthIXXGF\nfd9s/4XXc/f73T3f3fOzs7P3N06L9de3l5KWlMD1J+QFHUUk5owb2JnvHtuH/0xdw4S5G4KOI80k\n3EEwG4H3zOwN4PMbZ9z9rw2dZGbJ1BW/J9z9hdA5hfV+/gDwaujhOqD+bM25gD6JYVi4YQevfbaR\n68b1o2MbrfYgciBuOXkAM1Zt58fPz2NYTqZ6UuJAuC3AlcAkIIUwb4OwuiGIDwGL6hdKM+ta77Dz\ngPmh7ycAF5lZqpn1BvKAaWHmi2t3vrOUjLQkvn1Un6CjiMSs5NBN8slJCVz9hK4HxoOwWoDu/ssD\neO6xwDeBz8xsTmjfT4CLzWwEdd2bq4Dvhl5jgZk9AyykbgTpNe6uT2Aj5q0rZuLCQm4+qT+Z6clB\nxxGJad2yWvHXCw7hikdn8JvXFmrS7BausdUg7nD3Ow7kGHf/iH1f1/vSibPd/bfAbxt6Pflfd05c\nSlZ6MpeP7RV0FJEWYdzAzlx1TB/u/2AFR/XL5tShXYKOJBHSWAvw243c7mDU3bt3R5MlkrDNXL2d\nyUuK+NGpA8hIU+tPpKnccvIAPl2xlVufn8ew3ExysloFHUkioLFrgA/wv9f89t7ahI6RANw5cSkd\nWqdw6ZheQUcRaVFSkhL4+0WHUl1Ty41Pab7QlqrBFuABXvuTZjBj1TY+KtjC7WcMonVquIN5RSRc\nvTq25jfnDeWmp+fy93cLuPmk/kFHkiYW7ihQiTL3TC6gfesUvnFEz6CjiLRY5x2ay/kjc7j73WVM\nW7kt6DjSxFQAY9D89SVMXlLElUf1plWKpksViaRfnTOU7u3TuenpOZSUVwUdR5qQCmAMuu+95WSk\nJvHNMWr9iURam9Qk7rxwBJt27ObnL89v/ASJGeEuiNvfzCaZ2fzQ4+Fmdntko8m+LC/ayevzN/Kt\nI3vSViM/RZrFyB7tuH5cHi/P2cBLs9cHHUeaSLgtwAeAHwNVAO4+j7rbH6SZ3ffeclKTErhibO+g\no4jElWuO78thPdvxs5fms3ZbWdBxpAmEWwDT3X3vacmqmzqMNGzd9jJemr2ei0f1oIPm/BRpVkmJ\nCdx14Qgc+MEzc6mp3e+1ASTKhFsAt5hZX0KrM5jZV6mbIFua0QMfrMAMrjpGc36KBKF7+3TuOHsI\n01Zt48EPVwQdRw5SuAXwGuBfwEAzWw/cCHw/YqnkC7bvquTpGWs5d0QOXTM1K4VIUL4yModThnTm\nL28vZfGmBtcFlygXVgF09xXufiKQDQx096PcfVVEk8n/+M+0NeyuquU7av2JBMrM+N15w2jbKpmb\nnp5LRbXm7I9V4Y4C/Z2ZZbn7LncvNbN2ZvabSIeTOhXVNTz68SqO6Z9N/84NrkIlIs2gQ5tU/nD+\nMBZt3MFd7ywLOo4coHC7QE9z9+I9D9x9O3B6ZCLJ3l6Zu5Gi0gq+fZRGfopEixMHd+bC/O786/3l\nzFytWWJiUbgFMNHMPh92aGatAA1DbAbuzoMfrmBA5wyOzusYdBwRqednZw2mW1Yrbnl2HuWV6gqN\nNeEWwH8Dk8zsSjO7ApgIjI9cLNljSsFWFm8q5cqje2O2r+UVRSQobVKT+NNXh7Nyyy7+9NbioOPI\nfgp3EMyfqFuodhAwBPh1aJ9E2AMfrqBjm1TOGdEt6Cgisg9H9u3IpWN68siUVXy6YmvQcWQ/hD0X\nqLu/4e63uPsP3P2tSIaSOssKS3l/aRGXjulJapImvRaJVreeNpAe7dP50XPz2FWhOUJiRbijQM83\ns2VmVmJmO8ystJGV4qUJjP9kFSlJCXxjtCa9Folm6SlJ/Plrh7B2exl/eENdobEi3Bbgn4Cz3T3T\n3du6e4a7t41ksHhXuruKF2et56zh3WjfOiXoOCLSiFG923P5kb15/NPVfLJcXaGxINwCWOjuiyKa\nRP7HC7PWs6uyRkseicSQH54ygJ4d0rn1+XmUVaorNNqFWwBnmNnTZnZxqDv0fDM7P6LJ4pi78/in\nqxmem8mI7llBxxGRMLVKSeQP5w9nzbYy/vL20qDjSCPCLYBtgTLgZOCs0HZmpELFu0+Wb6Vg806+\nqWt/IjFnTN8OXDK6Bw9PWcnM1duDjiMNSArnIHe/PNJB5L8e/3Q1WenJnHWIbn0QiUW3nTaIyYuL\n+NFzc3nt+qNJS9Yo7mgU7ijQNDO7xszuNbOH92yRDhePNpaU8/bCQi7M765/NCIxqk1qEr87fxjL\ni3bxj3c1V2i0CrcL9HGgC3AK8D6QC5Q2dIKZdTezyWa2yMwWmNkNof3tzWxi6LaKiWbWLrTfzOzv\nZlZgZvPMbOSB/1qx68mpa6h15xtHqPtTJJYd2z+br4zM5V/vr2DhBt01Fo3CLYD93P1nwC53Hw+c\nAQxr5Jxq4AfuPggYDVxjZoOB24BJ7p4HTAo9BjgNyAttVwH37ddv0gJU1dTy5PS1HNc/mx4d0oOO\nIyIH6fYzBpGVnsxtL8zTCvJRKNwCWBX6WmxmQ4FMoFdDJ7j7RnefFfq+FFgE5ADn8N95RMcD54a+\nPwd4zOt8CmSZWddwf5GWYPLizRSVVvB1tf5EWoR2rVP4xVlDmLeuhEemrAw6juwl3AJ4f6ir8nZg\nArAQ+GO4L2JmvYBDgalAZ3ffCHVFEugUOiwHWFvvtHWhfXHj6elr6ZSRyvEDsoOOIiJN5MzhXTlh\nYCf+8vZS1m4rCzqO1BNuAZzk7tvd/QN37+PunYC3wznRzNoAzwM3untDHeH7WurgC30GZnaVmc0w\nsxlFRUVhhY8FhTt2M3nJZr5yWC5JiWFP0SoiUc7M+PW5Q0kw+MmLn+GurtBoEe7/tM/vY99zjZ1k\nZsmhc59w9xdCuwv3dG2Gvm4O7V8HdK93ei6wYe/ndPf73T3f3fOzs1tOS+m5meuodbggv3vjB4tI\nTOmW1YofnTqQD5dt4eU5X/hvTQLSYAE0s4Fm9hUgs/4MMGZ2GZDWyLkGPAQscve/1vvRBODS0PeX\nAi/X2/+t0GjQ0UDJnq7Slq621nlmxlqO6N2e3h1bBx1HRCLgktE9GdE9i1+/upDissqg4wiNtwAH\nUDfjSxb/nQHmLGAk8J1Gzh0LfBMYZ2ZzQtvpwB+Ak8xsGXBS6DHA68AKoAB4ALh6/3+d2DR15TZW\nby3jolFq/Ym0VIkJxu/OG0ZxeRW/f10rRkSDBmeCcfeXgZfNbIy7f7I/T+zuH7Hv63oAJ+zjeAeu\n2Z/XaCmenr6GjLQkThsaV4NeReLO4G5t+fZRvfnXBys4f2QOR/TpEHSkuBbuNcDzzKytmSWb2SQz\n22Jml0Q0WZwoKavijfmbOHdEjmZ+EYkDN5yYR267Vvzkxc+oqK4JOk5cC7cAnhwawXkmdYNV+gM/\njFiqODJh7noqqmu58HB1f4rEg/SUJH597lCWF+3iX++vCDpOXAu3ACaHvp4OPOnu2yKUJ+48P2s9\nA7tkMDQnM+goItJMjh/QiTOGdeXuyQWs3ror6DhxK9wC+IqZLQbygUlmlg3sjlys+LCiaCdz1hZz\n/si4ut9fRICfnTmY5ATj5y8v0L2BAQmrALr7bcAYIN/dq4Bd1E1dJgfhpTkbMINzRqgAisSbLplp\n3HzyAN5fWsSb8zcFHScuNTgK1MzGufu79Vd/r7u973MvfPEsCYe789Ls9Yzt25HObRu8pVJEWqhL\nx/Tk+Znr+OUrCzm6fzZtUsNaolWaSGMtwGNDX8/ax6YV4Q/CrDXbWbOtjPMOVetPJF4lJSbwm/OG\nUli6m7smLg06Ttxp7D7AX4S+akX4JvbCrPWkJSdwytAuQUcRkQCN7NGOiw7vwSMfr+Irh+UyqGvb\noCPFjca6QG9u6Od7TXEmYaqsruXVeRs5ZUgXdXmICLeeOoA352/k5y/P55nvjtn7UpNESGNdoBmh\nLR/4PnVPi/NjAAAb9ElEQVTLE+UA3wMGRzZayzV5yWZKyqs4V92fIgJkpadw66kDmb5qOy/NWR90\nnLjRYAF091+6+y+BjsBId/+Bu/8AOIy61RrkALw0ez0d26RwdL+OQUcRkShxQX53DumexW9fW8yO\n3VWNnyAHLdz7AHsA9acvr6SRFeFl30rKq5i0aDNnHdJN6/6JyOcSEoxfnzOErbsquGvisqDjxIVw\n/wd+HJhmZneY2S+oW9l9fORitVwTFxZSWVPL2Yd0CzqKiESZ4blZfH1UD8Z/sorFmxpaP1yaQrg3\nwv8WuBzYDhQDl7v77yMZrKV6bd4GcrJaMaJ7VtBRRCQK/fCUAbRNS9IMMc0g7D44d5/l7n8LbbMj\nGaqlKimr4sNlWzhzeFeN8hKRfcpKT+FHpw5k2sptvDIvLtYED4wuQjWjtxZuorrWOX2Y1v0TkS93\nQX53hua05XevLaKssjroOC2WCmAzev2zjeS2a8XwXK38ICJfLjHBuOOsIWzasZt7Jy8POk6LpQLY\nTIrLKvlo2RbOUPeniIQhv1d7zjs0h/s/WKElkyKkwQJoZqVmtmMfW6mZaYjSfnh7QSHVtc6ZwzT6\nU0TCc9tpA0lONH796qKgo7RIjd0In+HubfexZbi7JqzbD69+tpHu7VsxNEdvm4iEp3PbNK47IY93\nFhXy3pLNQcdpcfarC9TMOplZjz1bpEK1NNt3VTKlYAtnDOum7k8R2S+Xj+1Frw7p/PrVhVTV1AYd\np0UJqwCa2dlmtgxYCbwPrALeiGCuFuWtBZuoqXXOHK7RnyKyf1KTErn9jMEsL9rFvz9dHXScFiXc\nFuCvgdHAUnfvDZwATIlYqhbm9fmb6NE+nSHd1P0pIvvvhEGdODqvI3dOXMq2XZWNnyBhCbcAVrn7\nViDBzBLcfTIwIoK5Wowdu6v4ZPkWThnSWd2fInJAzIyfnTmYXZU13KmFc5tMuAWw2MzaAB8AT5jZ\n3wDdnRmG95cUUVXjnDxEC9+KyIHr3zmDS47owRNTV7NkU2nQcVqEcAvgOUAZcBPwJrAcOLOhE8zs\nYTPbbGbz6+27w8zWm9mc0HZ6vZ/92MwKzGyJmZ2y/79KdJq4sJAOrVMY2aNd0FFEJMbdeGJ/MtKS\n+dWrmie0KYRbAH/u7rXuXu3u493978CtjZzzKHDqPvbf6e4jQtvrAGY2GLgIGBI6514zSwwzW9Sq\nqqll8pLNjBvYicQEdX+KyMFp1zqFm07MY0rBViYt0m0RByvcAnjSPvad1tAJ7v4BsC3M5z8HeMrd\nK9x9JVAAjArz3Kg1dcU2SndXq/tTRJrMN0b3pG92a373+iLdFnGQGpsJ5vtm9hkwwMzm1dtWAvMO\n8DWvDT3Hw2a2p18wB1hb75h1oX0x7e2Fm0hLTuAorfwuIk0kOTGBn54xiBVbdFvEwWqsBfgf4Cxg\nQujrnu0wd7/kAF7vPqAvdSNINwJ/Ce3fV//gPju4zewqM5thZjOKiooOIELzcHfeWVjI0XnZtEqJ\n+d5cEYkixw/oxFH9OnLXO8soLtNtEQeqsanQStx9lbtfTF2rrIq6wtTmQGaCcfdCd69x91rgAf7b\nzbkO6F7v0Fxgw5c8x/3unu/u+dnZ2fsbodks2LCDDSW7OWlw56CjiEgLY2bcfuYgSndX8fdJBUHH\niVnhzgRzLVAITAReC22v7u+LmVn9qVDOA/aMEJ0AXGRmqWbWG8gDpu3v80eTtxcWkmBwwsBOQUcR\nkRZoYJe2XHh4dx77ZBUrinYGHScmhTsI5kZggLsPcfdhoW14QyeY2ZPAJ9RdP1xnZlcCfzKzz8xs\nHnA8dbdV4O4LgGeAhdTdZnGNu9cc4O8UFd5esIn8nu3p0CY16Cgi0kLdfNIAUpMS+P0bi4OOEpOS\nwjxuLVCyP08c6jbd20MNHP9b4Lf78xrRau22MhZvKuWnpw8KOoqItGDZGalcfXw//u+tJUxdsZUj\n+nQIOlJMCbcFuAJ4L3Sz+s17tkgGi2WTQ8uWnKjrfyISYVeM7U3XzDR+9/oiamt1c/z+CLcArqHu\n+l8KkFFvk32YvHgzPTuk07tj66CjiEgL1yolkVtOHsDcdSW8Mm+fYwflS4TVBeruvwQws4y6h64r\nrl9id1UNn6zYyoX53Rs/WESkCZx3aA4PfbSSP725hFOGdCEtWbdehSPcUaBDzWw2daM2F5jZTDMb\nEtlosWnqym3srqrluAEa/SkizSMhwbj9jEGsLy5n/Mergo4TM8LtAr0fuNnde7p7T+AH1N3HJ3t5\nb8lmUpISGK2L0SLSjI7s15FxAztx9+QCtmvNwLCEWwBbh9YABMDd3wN0gWsf3l9SxJg+HTT7i4g0\nux+fNpBdFdXcPVk3x4cj7FGgZvYzM+sV2m4HVkYyWCxas7WMFVt2cdyA6J2hRkRarrzOGXztsLqb\n49duKws6TtQLtwBeAWQDL4S2jsDlkQoVq95bWnf7g67/iUhQbjqpP4kJxv+9tSToKFEv3AI4FLjJ\n3UeGthuB3hHMFZPeW1Kk2x9EJFBdMtO48qjeTJi7gc/W7df8JXEn3AL4FvCumdW/s/vBCOSJWbur\navh4+RaO66/uTxEJ1neP7Uv71in8/o1FWjm+AeEWwCXA/1E3G8yRoX1a4ryeabr9QUSiRNu0ZK4b\n14+Pl2/l/aXRu2xc0MItgO7urwJnA3eHVofQnxX1vLekSLc/iEjU+MYRPenRPp0/vLFYU6R9iXAL\noAG4+zLgaOAYoMHVIOLN+0s3M1q3P4hIlEhJSuCWUwaweFMpL89dH3ScqBRWAXT3Q+t9v8vdLwD6\nRCxVjNlUspvlRbs4ul/HoKOIiHzuzGFdGdKtLX9+aykV1TG9wlxENFgAzexHoa//MLO/19+AW5ol\nYQyYUrAFgCP7qftTRKJHQoJx66kDWV9czhOfrgk6TtRpbDLsRaGvMyIdJJZNWb6F9q1TGNSlbdBR\nRET+x9F5HTmybwfunlzA1/JzyUhLDjpS1GiwBejur5hZIjDU3cfvvTVTxqjm7nxcsJUxfTqQkKCB\nsSISXczqWoHbdlXywIeawKu+Rq8BunsNcFgzZIlJK7bsYtOO3Yzpq+5PEYlOh3TP4oxhXXnwwxUU\nlVYEHSdqhDsKdLaZTTCzb5rZ+Xu2iCaLER+Hrv+N1QAYEYlit5wygIrqWu7RRNmfC7cAtge2AuOA\ns0LbmZEKFUumFGylW2YavTqkBx1FRORL9e7Ymgvyc3li6mpNlB0S7orwmvh6H2prnU9WbOWkwZ0x\n0/U/EYlu15+Qx/Oz1nPXO8v4ywWHBB0ncOGuCJ9rZi+a2WYzKzSz580sN9Lhot3CjTsoKa9irG5/\nEJEY0DWzFZeO6ckLs9extLA06DiBC7cL9BFgAtANyAFeCe2La5/f/9dX1/9EJDZ8/7h+tE5J4s9a\nLinsApjt7o+4e3Voe5S69QHj2pTlW+nXqQ2d26YFHUVEJCztW6fwnaP78PbCQmav2R50nECFWwC3\nmNklZpYY2i6hblBM3KqsrmX6ym2M1e0PIhJjrjy6Nx1ap8T9orn7syL8BcAmYCPw1dC+uDV7zXbK\nq2o4Urc/iEiMaZOaxNXH1y2XtOdSTjwKdzLsNe5+trtnu3sndz/X3Vc3dI6ZPRwaNDO/3r72ZjbR\nzJaFvrYL7bfQHKMFZjbPzEYe3K8VeZ+s2IoZjO6tFqCIxJ5vHNGDrplp/N9bS+J20dzGJsP+wiTY\ne02I3ZBHgVP32ncbMMnd84BJoccApwF5oe0q4L79/UWa2/RV2xjYpS2Z6ZpXT0RiT1pyIteNy2PO\n2mLeXbw56DiBaKwFOAOYGdrOrvf9nu1LufsHwLa9dp8D7JlDdDxwbr39j3mdT4EsM+sa7i/R3Kpq\napm1uphRvdoFHUVE5IB9LT+Xnh3S+fPbS+Ny0dwGb4SvP+G1md3YBBNgd3b3jaHn3mhmnUL7c4C1\n9Y5bF9q3ce8nMLOrqGsl0qNHj4OMc2AWbNhBeVUNh/duH8jri4g0heTEBG48MY+bnp7L6/M3cubw\nbkFHalbhDoIBiOSfB/uaRmWfr+fu97t7vrvnZ2cHcyfG9JV1DdtRvVQARSS2nX1IDnmd2vDXiUup\nrqkNOk6z2p8C2BQK93Rthr7u6XheB3Svd1wusKGZs4Vt2qpt9OyQTifd/yciMS4xwfjByQNYUbSL\nF2evDzpOs2psEEypme0wsx3A8D3f79l/AK83Abg09P2lwMv19n8rNBp0NFCyp6s02tTWOjNWbeNw\ntf5EpIU4ZUhnhuVk8rdJy6isjp9WYGML4ma4e9vQllTv+wx3b3D5czN7EvgEGGBm68zsSuAPwElm\ntgw4KfQY4HVgBVAAPABcfZC/V8QsL9rJ9rIqdX+KSIthZtx8cn/WbS/n2ZlrGz+hhQhrNYgD4e4X\nf8mPTtjHsQ5cE6ksTWnaqrrrfxoAIyItyXH9sxnZI4u73y3gKyNzSUtODDpSxDX3NcCYN33lNjq2\nSdX6fyLSopjVXQvcWLKbJ6etCTpOs1AB3E/TV21nVO92Wv9PRFqcI/t2YHSf9twzeTnllTVBx4k4\nFcD9sL64nPXF5RoAIyIt0p5W4JadFTz+6aqg40ScCuB+2HP/nwqgiLRUh/dqzzH9s7nvveXsrKgO\nOk5EqQDuh2mrtpGRmsSgrg0OgBURiWk3n9Sf7WVVjP94VdBRIkoFcD9MX7mNkT3bkZig638i0nKN\n6J7FuIGduP+DFZTurgo6TsSoAIapuKySZZt3Mkq3P4hIHLjpxP6UlFfx6JRVQUeJGBXAMM1ZWwzA\nod2zAk4iIhJ5w3IzOWlwZx74cAUl5S2zFagCGKY5a4sxg+EqgCISJ248MY8du6t5+KOVQUeJCBXA\nMM1ZW0z/Thm0SY3Y5DkiIlFlSLdMTh3ShYc/WklJWctrBaoAhsHdmbu2mBFq/YlInLnhxDxKK6p5\n8KMVQUdpciqAYVi9tYztZVWM6KECKCLxZVDXtpwxrCuPTFlFcVll0HGalApgGGav3Q6gFqCIxKXr\nT8hjZ0U1D7Wwa4EqgGGYs6aY9JRE+nfOCDqKiEizG9Alo0W2AlUAwzBnbTHDczN1A7yIxK2W2ApU\nAWzE7qoaFm7cwYju7YKOIiISmJbYClQBbMTCjTuoqnFd/xORuNfSWoEqgI2YvSY0A4xGgIpInGtp\nrUAVwEbMWVtM18w0OrdNCzqKiEjgWlIrUAWwEXPWblfrT0QkZECXDE4b2oVHp6yK+dlhVAAbsHVn\nBWu3lev6n4hIPdefUDc7zENTYrsVqALYgD0rQGgEqIjIfw3q2pZThnTmkSkrY3qlCBXABsxeU0xi\ngjEsJzPoKCIiUeX6E/Io3V0d0+sFqgA24LP1JeR1akOrlMSgo4iIRJUh3TI5cVBnHvpoBTtidNV4\nFcAv4e7MX1+i1p+IyJe44YS69QLHx2grMJACaGarzOwzM5tjZjNC+9qb2UQzWxb6GuiFt8IdFWzd\nVcmQbm2DjCEiErWG5WZywsBOPPjRSkpjsBUYZAvweHcf4e75oce3AZPcPQ+YFHocmAUbSgAYqhag\niMiXuu6EPErKq3j809VBR9lv0dQFeg4wPvT9eODcALMwf/0OzOpGO4mIyL6N6J7Fsf2zefDDlZRV\nVgcdZ78EVQAdeNvMZprZVaF9nd19I0Doa6d9nWhmV5nZDDObUVRUFLGACzaU0Ltja1qnJkXsNURE\nWoLrT+jHtl2VPPHpmqCj7JegCuBYdx8JnAZcY2bHhHuiu9/v7vnunp+dnR2xgAs27GBoN3V/iog0\n5rCe7RnbrwP/+mAFu6tqgo4TtkAKoLtvCH3dDLwIjAIKzawrQOjr5iCyAWzfVcn64nINgBERCdN1\n4/LYsrOCJ6fFTiuw2QugmbU2s4w93wMnA/OBCcClocMuBV5u7mx7LNiwA9AAGBGRcI3u04FRvdvz\nz/eXx0wrMIgWYGfgIzObC0wDXnP3N4E/ACeZ2TLgpNDjQMwPjQBVC1BEJHw3nJBH4Y4Knp25Lugo\nYWn2ER7uvgI4ZB/7twInNHeefVmwYQc5Wa3ISk8JOoqISMw4sm8HRvbI4p/vLefC/O6kJEXTjQZf\nFN3pArJgfQlDc9T6ExHZH2bGdePyWF9czouzo78VqAK4l50V1azcuoshGgEqIrLfjhuQzbCcTO59\nbznVNbVBx2mQCuBeFm3cgTtqAYqIHAAz49px/Vi9tYxX5m0IOk6DVAD3Mn99aAo0tQBFRA7ISYM6\nM6BzBne/W0BNrQcd50upAO5lwYYddGyTSqe2aUFHERGJSQkJda3A5UW7eGP+xqDjfCkVwL3M1wAY\nEZGDdvqwrvTJbs3d7xZQG6WtQBXAenZX1VCweafu/xMROUiJCcY1x/Vj8aZSJi0ObGKvBqkA1rOs\ncCfVta4RoCIiTeCcEd3o3r4Vd7+7DPfoawWqANazeFPdFGhaAklE5OAlJSbw/WP7MXddCR8u2xJ0\nnC9QAaxnaWEpqUkJ9GifHnQUEZEW4SuH5dClbRp3Ty4IOsoXqADWs6RwJ3md25CYYEFHERFpEVKT\nEvnusX2YtnIbU1dsDTrO/1ABrGfpplL6d84IOoaISIty0eE96NgmJepagSqAISVlVWzasZsBKoAi\nIk2qVUoi3z66Dx8u28LctcVBx/mcCmDI0s2lAPTvogIoItLULhndk8xWyfzj3ehpBaoAhizZVFcA\n1QIUEWl6bVKTuHxsL95ZVPj5iPugqQCGLC0sJSM1ia6ZmgJNRCQSLjuyF61TErln8vKgowAqgJ9b\nvKmU/l0yMNMIUBGRSMhKT+GSMT15bd4GVm7ZFXQcFUAAd2dpYSkDdP1PRCSivn1UH5ITE7jvveCv\nBaoAAkWlFRSXVen6n4hIhGVnpHLxqB68MGs967aXBZpFBRBYUhgaAaoCKCIScVcd0wczuP+DFYHm\nUAHkvyNA+3duE3ASEZGWr1tWK84/NJenpq9lc+nuwHKoAFI3ArRjm1Q6tEkNOoqISFz4/nF9qa6p\n5aEPVwaWQQWQujlAB3RR609EpLn06tiasw7pxr8/XU1xWWUgGeK+ANbWOssKNQeoiEhzu/q4fuyq\nrOGRKasCef24L4Dri8spq6zRCFARkWY2oEsGJw/uzKMfr2JnRXWzv37UFUAzO9XMlphZgZndFunX\nW7xJc4CKiATl2nH9KCmv4t+frm72146qAmhmicA9wGnAYOBiMxscyddcGroFIq+TrgGKiDS34blZ\nHJ3XkQc/XMnuqppmfe2oKoDAKKDA3Ve4eyXwFHBOJF9wyaZScrJakZGWHMmXERGRL3Ht8f3YsrOC\np6evbdbXjbYCmAPUfwfWhfZ9zsyuMrMZZjajqKjooF9QU6CJiATriD4dOLxXO8Z/vAp3b7bXTWq2\nVwrPvmai/p93w93vB+4HyM/PP+h36tHLR1HezM1uERH5X787bxiZ6cnNuiBBtBXAdUD3eo9zgQ2R\nfMEuWv5IRCRweQGMxI+2LtDpQJ6Z9TazFOAiYELAmUREpAWKqhagu1eb2bXAW0Ai8LC7Lwg4loiI\ntEBRVQAB3P114PWgc4iISMsWbV2gIiIizUIFUERE4pIKoIiIxCUVQBERiUsqgCIiEpdUAEVEJC6p\nAIqISFyy5px4tKmZWRHQFItIdQS2NMHzxCO9dwdO792B03t34OLhvevp7tmNHRTTBbCpmNkMd88P\nOkcs0nt34PTeHTi9dwdO791/qQtURETikgqgiIjEJRXAOvcHHSCG6b07cHrvDpzeuwOn9y5E1wBF\nRCQuqQUoIiJxSQVQRETiUlwXQDM71cyWmFmBmd0WdJ5oZmbdzWyymS0yswVmdkNof3szm2hmy0Jf\n2wWdNVqZWaKZzTazV0OPe5vZ1NB797SZpQSdMRqZWZaZPWdmi0OfvzH63IXHzG4K/Xudb2ZPmlma\nPnf/FbcF0MwSgXuA04DBwMVmNjjYVFGtGviBuw8CRgPXhN6v24BJ7p4HTAo9ln27AVhU7/EfgTtD\n79124MpAUkW/vwFvuvtA4BDq3kN97hphZjnA9UC+uw8FEoGL0Ofuc3FbAIFRQIG7r3D3SuAp4JyA\nM0Utd9/o7rNC35dS959QDnXv2fjQYeOBc4NJGN3MLBc4A3gw9NiAccBzoUP03u2DmbUFjgEeAnD3\nSncvRp+7cCUBrcwsCUgHNqLP3efiuQDmAGvrPV4X2ieNMLNewKHAVKCzu2+EuiIJdAouWVS7C/gR\nUBt63AEodvfq0GN9/vatD1AEPBLqPn7QzFqjz12j3H098GdgDXWFrwSYiT53n4vnAmj72Kd7Qhph\nZm2A54Eb3X1H0HligZmdCWx295n1d+/jUH3+vigJGAnc5+6HArtQd2dYQtdFzwF6A92A1tRd8tlb\n3H7u4rkArgO613ucC2wIKEtMMLNk6orfE+7+Qmh3oZl1Df28K7A5qHxRbCxwtpmtoq6rfRx1LcKs\nUNcU6PP3ZdYB69x9aujxc9QVRH3uGncisNLdi9y9CngBOBJ97j4XzwVwOpAXGhGVQt3F4QkBZ4pa\noWtWDwGL3P2v9X40Abg09P2lwMvNnS3aufuP3T3X3XtR9zl7192/AUwGvho6TO/dPrj7JmCtmQ0I\n7ToBWIg+d+FYA4w2s/TQv989750+dyFxPROMmZ1O3V/iicDD7v7bgCNFLTM7CvgQ+Iz/Xsf6CXXX\nAZ8BelD3D+5r7r4tkJAxwMyOA25x9zPNrA91LcL2wGzgEnevCDJfNDKzEdQNHkoBVgCXU/fHuz53\njTCzXwIXUjeKezbwbequ+elzR5wXQBERiV/x3AUqIiJxTAVQRETikgqgiIjEJRVAERGJSyqAIiIS\nl1QARUQkLqkAiohIXFIBFIkCobUC/xZau+2z0E3y9X/ey8zKzWzOl5x/h5ndcoCv3crM5phZpZl1\nPJDnEIlFKoAi0eHHwAp3HwL8Hbh6H8csd/cRTf3C7l4eet64nRNS4pMKoEjAQsv7nOfufwvtWgn0\nC+O8n5rZEjN7BxhQb/8lZjYt1Kr7V2jxZ8zsZ6FV1SeGVgc/oBajSEuR1PghIhJhJwLd63Vvtgfe\naegEMzuMuom1D6Xu3/EsYKaZDaJu7sex7l5lZvcC3zCzhcBX9j4+Er+MSKxQARQJ3gjg5+7+TwAz\nexCY18g5RwMvuntZ6Jw9K5mcABwGTK9bAIBW1C0V1B542d3LQ8e/0tS/hEisUReoSPDaAXsKWRJw\nMhBOgdrXTPYGjHf3EaFtgLvfwb4X4BWJayqAIsFbCowOfX8T8Jq7r2zknA+A80IjODOAs0L7JwFf\nNbNOAGbW3sx6Ah8BZ5lZmpm1Ac5o8t9CJMaoC1QkeE8Cb5hZAfAJcFVjJ7j7LDN7GpgDrKZurUbc\nfaGZ3Q68bWYJQBVwjbt/GuomnRs6fgZQEpHfRiRGaD1AkRhgZr2AV9196EE8Rxt332lm6dS1IK9y\n91n1fr4KyHf3LQcZVyQmqAtUJDbUAJlfdiN8mO4PnT8LeH5P8dtzIzyQDNQefFSR2KAWoIiIxCW1\nAEVEJC6pAIqISFxSARQRkbikAigiInFJBVBEROKSCqCIiMQlFUAREYlL/w+S8xmkl+BOUwAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d4be810>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Now use the method to make a plot of horizontal distance versus theta; again use default values\n",
    "theta = np.linspace(0, np.pi/2.0,100)\n",
    "xh = [hdist(th) for th in theta]\n",
    "plt.plot(theta*180/np.pi, xh)\n",
    "plt.xlabel(r'$\\theta$ [deg]')\n",
    "plt.ylabel('Horizontal distance [m]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Understanding When Solutions Exist\n",
    "\n",
    "We have now completed step 1: we have a function that returns the horizontal distance travelled as a function of angle. The next step is to write a function that will use this to solve for the angle as a function of target distance. One important thing that the plot above should make clear is that there will in general be either two solutions or zero solutions. There is some maximum distance which the arrow can travel, which in this example is just over 350 m, for an arow fired at just under $40^\\circ$. If the target distance is smaller than this there will be two solutions, one corresponding to a high arc and the other to a shallow arc. If the target distance is larger, the arrow cannot reach it at all. This illustrates an important point about ODE systems of this type: there are no guarantees that solutions exist, or that they are unique. This is why physical insight is crucial.\n",
    "\n",
    "Having seen this, we can now go about designing a shooting method that will solve the problem for us. This will take part in two steps. First, find the maximum distance for a given set of input parameters. Second, if solutions exist, search for the two solutions.\n",
    "\n",
    "For the first part of this, we will use the scipy library function minimize_scalar, which finds minima (or maxima if we just apply a minus sign) of a scalar function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import the minimize_scalar function\n",
    "from scipy.optimize import minimize_scalar\n",
    "\n",
    "# Define a function that is just the negative of our existing one, so we can minimize it\n",
    "def hdist_neg(theta, v=100.0, dt=10.0, dt_lim=1.0e-10, z_lim=1.0e-6, verbose=False):\n",
    "    return -hdist(theta, v=v, dt=dt, dt_lim=dt_lim, z_lim=z_lim, verbose=verbose)\n",
    "\n",
    "# Define a function that returns the angle of the maximum and the maximum distance\n",
    "def hdist_max(v=100.0, dt=10.0, dt_lim=1.0e-10, z_lim=1.0e-6, verbose=False):\n",
    "    \n",
    "    # Find minimum; we can speed this up using the bracket option, which lets us provide\n",
    "    # a search region by specifying values (a, b, c) such that hdist(b) < hdist(a) and\n",
    "    # hdist(b) < hdist(c)\n",
    "    result = minimize_scalar(hdist_neg, bracket=(0, np.pi/4, np.pi/2))\n",
    "    \n",
    "    # Return result\n",
    "    theta = result.x\n",
    "    dist = -result.fun\n",
    "    return theta, dist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "theta_max = 37.223351 deg\n",
      "maximum distance = 360.971862 m\n"
     ]
    }
   ],
   "source": [
    "# Demonstrate the method for the default arguments\n",
    "theta_max, hd_max = hdist_max()\n",
    "print(\"theta_max = {:f} deg\".format(theta_max*180/np.pi))\n",
    "print(\"maximum distance = {:f} m\".format(hd_max))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Putting it Together: the Shooting Method\n",
    "\n",
    "Now we are in a position to perform the final step: writing a function that will find the angle at which the crossbow should be fired to hit a target at a specified distance. The function will check if the range given is reachable for any angle. If it is, it will return the two possible angles that can hit that distance. If not, it will return an empty solution. We will again use a scipy library function, brentq."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import brentq\n",
    "from scipy.optimize import brentq\n",
    "\n",
    "# Define a function for the difference between the distance of the shot and the target\n",
    "def err(theta, target_dist, v, dt, dt_lim, z_lim, verbose):\n",
    "    shot_dist = hdist(theta, v=v, dt=dt, dt_lim=dt_lim, z_lim=z_lim, verbose=False)\n",
    "    if verbose:\n",
    "        print(\"theta = {:f}, shot distance = {:f}, miss distance = {:e}\"\n",
    "              .format(theta, shot_dist, shot_dist-target_dist))\n",
    "    return shot_dist - target_dist\n",
    "\n",
    "# Define the function\n",
    "def target_angle(dist, v=100.0, dt=10.0, dt_lim=1.0e-10, z_lim=1.0e-6, verbose=False):\n",
    "    \n",
    "    # First find the maximum distance\n",
    "    theta_max, dist_max = hdist_max(v=v, dt=dt, dt_lim=dt_lim, z_lim=z_lim)\n",
    "    \n",
    "    # Check if a solution exists; if not, return empty value\n",
    "    if dist > dist_max:\n",
    "        return None\n",
    "    \n",
    "    # If we are here, two solutions exist, one with theta < theta_max, one with\n",
    "    # theta > theta_max; find both\n",
    "    if verbose:\n",
    "        print(\"Finding low angle:\")\n",
    "    theta1 = brentq(err, 0.0, theta_max,\n",
    "                    args=(dist, v, dt, dt_lim, z_lim, verbose))\n",
    "    if verbose:\n",
    "        print(\"\")\n",
    "        print(\"Finding high angle:\")\n",
    "    theta2 = brentq(err, theta_max, np.pi/2, \n",
    "                    args=(dist, v, dt, dt_lim, z_lim, verbose))\n",
    "\n",
    "    # Return result\n",
    "    return theta1, theta2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finding low angle:\n",
      "theta = 0.000000, shot distance = 0.030516, miss distance = -2.999695e+02\n",
      "theta = 0.649670, shot distance = 360.971862, miss distance = 6.097186e+01\n",
      "theta = 0.539925, shot distance = 354.757392, miss distance = 5.475739e+01\n",
      "theta = 0.269962, shot distance = 275.407603, miss distance = -2.459240e+01\n",
      "theta = 0.353630, shot distance = 311.730895, miss distance = 1.173089e+01\n",
      "theta = 0.326609, shot distance = 301.373569, miss distance = 1.373569e+00\n",
      "theta = 0.323215, shot distance = 299.985205, miss distance = -1.479479e-02\n",
      "theta = 0.323251, shot distance = 300.000106, miss distance = 1.059520e-04\n",
      "theta = 0.323251, shot distance = 299.999998, miss distance = -1.604579e-06\n",
      "theta = 0.323251, shot distance = 300.000001, miss distance = 8.741292e-07\n",
      "theta = 0.323251, shot distance = 299.999998, miss distance = -1.803546e-06\n",
      "theta = 0.323251, shot distance = 300.000001, miss distance = 9.095139e-07\n",
      "theta = 0.323251, shot distance = 300.000001, miss distance = 9.460168e-07\n",
      "theta = 0.323251, shot distance = 299.999998, miss distance = -1.821797e-06\n",
      "theta = 0.323251, shot distance = 299.999998, miss distance = -1.833811e-06\n",
      "theta = 0.323251, shot distance = 299.999998, miss distance = -1.837926e-06\n",
      "theta = 0.323251, shot distance = 300.000001, miss distance = 9.467385e-07\n",
      "theta = 0.323251, shot distance = 299.999998, miss distance = -1.838627e-06\n",
      "theta = 0.323251, shot distance = 300.000001, miss distance = 9.469765e-07\n",
      "theta = 0.323251, shot distance = 300.000001, miss distance = 9.472082e-07\n",
      "theta = 0.323251, shot distance = 300.000001, miss distance = 9.473239e-07\n",
      "\n",
      "Finding high angle:\n",
      "theta = 0.649670, shot distance = 360.971862, miss distance = 6.097186e+01\n",
      "theta = 1.570796, shot distance = 0.000000, miss distance = -3.000000e+02\n",
      "theta = 0.805258, shot distance = 349.356157, miss distance = 4.935616e+01\n",
      "theta = 1.188027, shot distance = 229.977096, miss distance = -7.002290e+01\n",
      "theta = 0.963510, shot distance = 315.056927, miss distance = 1.505693e+01\n",
      "theta = 1.020686, shot distance = 297.359631, miss distance = -2.640369e+00\n",
      "theta = 1.012156, shot distance = 300.175048, miss distance = 1.750481e-01\n",
      "theta = 1.012686, shot distance = 300.001782, miss distance = 1.781857e-03\n",
      "theta = 1.012692, shot distance = 300.000000, miss distance = 1.915009e-08\n",
      "theta = 1.012692, shot distance = 300.000000, miss distance = -4.799688e-09\n",
      "theta = 1.012692, shot distance = 300.000000, miss distance = 5.684342e-14\n",
      "theta = 1.012692, shot distance = 300.000000, miss distance = -4.091589e-10\n",
      "Firing angles = 18.520918, 58.022963 deg\n"
     ]
    }
   ],
   "source": [
    "# Demonstrate the result\n",
    "th1, th2 = target_angle(300, verbose=True)\n",
    "print(\"Firing angles = {:f}, {:f} deg\".format(th1*180/np.pi, th2*180/np.pi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This method can clearly be generalised to any sort of complex ODE, just by changing the function \"err\" to something appropriate. For example, if we wanted to constrain the total path length travelled by the arrow, we could replace the function in our program by a function that integrated the path length returned by the ODE solver, and use the same basic structure. However, as always, remember that in general there are no guarantees that solutions exist, or how many there might be. This depends on the nature of the problem, and often requires some sort of physical insight into what is going on.\n",
    "\n",
    "One more note is that there is a large class of problems where, rather than shooting from one side and trying to hit the constraint, it is preferable to shoot from both sides and try to hit one or more points in the middle. For example, rather than integrating from the crossbow to the target, we could try to integrate the trajectory both from the crossbow forward and from the target backward, and iterate until the two solutions match in the middle. This is a bad idea for this particular system, but is a good idea in many cases, particularly when the \"target\" one is trying to hit represents a singular point of the solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finite Difference / Relaxation Methods\n",
    "\n",
    "Shooting should usually be your first go-to for most complex ODEs, but there are situations where other methods are preferable. One particularly important application is cases where we want to generate solutions many times for a range of parameter values. In our archery example, we might be concerned with a case where the goal is to hit 1000 targets, all placed at slightly different distances. Shooting is sub-optimal for this case because it does a lot of redundant work. Each time we shoot we have to integrate the ODEs describing the trajectory from the firing point to the point where the arrow impacts, so it costs the same computational expense to hit each target. However, this is inefficient. If two targets are 1 cm apart, then the trajectories the arrows are going to take to them are nearly identical. We should be able to take advantage of that information somehow.\n",
    "\n",
    "One way of handling a situation like this is known as either the relaxation method or the finite difference method. The idea of this method is that we take the ODE system and discretize it onto a mesh, and linearize if necessary, so that the ODE system plus the boundary conditions together form a system of linear algebraic equations. Such systems can be represented as linear matrix equations, and there are numerous methods of solving such systems. We use whichever ones we like. Often it is preferable not to solve the matrix equation exactly (since exact solutions of big matrix equations are expensive), but instead to approximate the solution iteratively. Iteration is also required if we have had to linearize a set of non-linear algebraic equations. We start with an initial guess at the solution, and then solve that linear system iteratively until we have something close enough to the right solution. This is why the method is known as relaxation: we start with a guess and then gradually relax it to the right solution.\n",
    "\n",
    "The main advantage of this is that, if we start with a guess that is close to correct, the iteration process will converge very quickly, requiring little work. If we are solving ODEs for a range of parameters, and the parameters are changing slowly, then we have a natural source of good guesses: our guess for a particular parameter is just our solution for the previous one. In situations like this, the first solution is often expensive to get, but the subsequent ones are dirt cheap, so for a large number of parameters we can win by a lot over shooting. A secondary advantage is that relaxation methods are significantly easier and more efficient to parallelise.\n",
    "\n",
    "For this reason relaxation methods are used extensively in astrophysics, in stellar evolution for example; the central part of the [MESA stellar evolution code](http://mesa.sourceforge.net/), for example, is a method of this type. The old state of the star is a good guess for the new state, so for stellar evolution calculations relaxation is much more numerically efficient than shooting."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "All these steps are easiest to understand via example; this example is somewhat artificial, but it nicely illustrates the method and its advantages. Let us consider a rod of heat-conducting metal. The two ends of the rod are placed in thermal baths that keep them at some fixed temperature, and heat conducts through the rod, so the temperature varies across it. We want to solve for the temperature of the rod versus position, and we want to find this for a large number of different temperature differences across the rod.\n",
    "\n",
    "Mathematically, this system obeys the time-independent diffusion equation, which in this 1D case can be cast as a 2nd order ODE:\n",
    "\\begin{equation}\n",
    "0 = \\frac{d}{dx}\\left[K(T) \\frac{dT}{dx}\\right] = K(T) \\frac{d^2 T}{dx^2} + K'(T) \\left(\\frac{dT}{dx}\\right)^2\n",
    "\\end{equation}\n",
    "where $K(T)$ is the thermal diffusivity, which in general can be a function of $T$; we assume that $K(T)$ is known. If we adopt units whereby the rod has length 1 and the temperature at the left hand side is 1, then the boundary conditions are $T = 1$ at $x=0$ and $T = T_f$ at $x=1$. If $K$ is a constant this equation is trivial: the solution is just that $T$ is a linear function going from $1$ to $T_f$ as $x$ goes from 0 to 1. However, if $K$ is not constant than the equation is non-linear and a numerical solution will generally be required. Our goal is to generate the solution for a wide range of $T_f$ values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Discretization\n",
    "\n",
    "To solve a problem like our example with a finite difference method, the first step is to discretize the equations to produce a system of algebraic equations from the ODEs. Consider a set of $N$ points distributed uniformly from $x=0$ to $x=1$, so point $i$ is located at position $x_i = i/(N-1)$ for $i=0 \\ldots N-1$. The spacing between the points is $\\Delta x = 1/N$. The next step is to write out approximations for the derivatives on evaluated at the mesh points. One might be tempted just to take\n",
    "\\begin{equation}\n",
    "\\left(\\frac{dT}{dx}\\right)_i \\approx \\frac{x_{i+1} - x_i}{\\Delta x},\n",
    "\\end{equation}\n",
    "but there is a better choice to be made. To see why, let us use our trick of writing down Taylor expansions again. The expansion of $T(x)$ about $x=x_i$ is\n",
    "\\begin{equation}\n",
    "T(x) = T_i + (x - x_i) T'_i + \\frac{(x - x_i)^2}{2} T''_i + O\\left((x-x_i)^3\\right),\n",
    "\\end{equation}\n",
    "where for compactness we have used $T_i$ as shorthand for $T(x_i)$ and similarly for the derivative terms. Re-arranging,\n",
    "\\begin{equation}\n",
    "T'_i = \\frac{1}{x-x_i}\\left[T(x) - T_i - \\frac{(x-x_i)^2}{2} T''_i + O\\left((x-x_i)^3\\right)\\right]\n",
    "\\end{equation}\n",
    "Now let us evaluate this for both $x=x_{i+1}$ and $x=x_{i-1}$:\n",
    "\\begin{eqnarray}\n",
    "T'_i & = & \\frac{1}{\\Delta x}\\left[T_{i+1} - T_i - \\frac{\\Delta x^2}{2} T''_i + \n",
    "O\\left(\\Delta x^3\\right)\\right] \\\\\n",
    "T'_i & = & -\\frac{1}{\\Delta x}\\left[T_{i-1} - T_i - \\frac{\\Delta x^2}{2} T''_i +\n",
    "O\\left(\\Delta x^3\\right)\\right]\n",
    "\\end{eqnarray}\n",
    "Adding these two equations, we have\n",
    "\\begin{equation}\n",
    "T'_i = \\frac{T_{i+1} - T_{i-1}}{2\\Delta x} + O\\left(\\Delta x^2\\right).\n",
    "\\end{equation}\n",
    "Thanks to the exact cancelation of the $T''_i$ terms, this approximation is accurate to order $\\Delta x^2$, while the more naive one above is only accurate to order $\\Delta x$. Analogous expansion of the second derivative gives the second-order accurate approxiation\n",
    "\\begin{equation}\n",
    "T''_i = \\frac{T_{i+1} - 2 T_i + T_{i-1}}{\\Delta x^2} + O\\left(\\Delta x^2\\right).\n",
    "\\end{equation}\n",
    "\n",
    "Armed with this approximation, we can rewrite our ODE as\n",
    "\\begin{equation}\n",
    "\\frac{T_{i+1} - 2 T_i + T_{i-1}}{\\Delta x^2} + \\frac{K'_i}{K_i} \\left(\\frac{T_{i+1}-T_{i-1}}{2\\Delta x}\\right)^2 = 0\n",
    "\\end{equation}\n",
    "where we have used the shorthand $K_i = K(T_i)$ and similarly for the $K'_i$ term.\n",
    "\n",
    "The boundary conditions can then be expressed naturally as\n",
    "\\begin{eqnarray}\n",
    "T_0 & = & 1 \\\\\n",
    "T_1 & = & T_f.\n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linearization\n",
    "\n",
    "The equations we have just derived are discrete, but they are not yet linear, since $T_{i+1}$ and $T_{i-1}$ appear squared in the second term, and the diffusivity $K$ and its derivatives also involve arbitrary functions of $T_i$. Our next step is to find a way of approximating this by a linear system. Just as there is more than one way to discretize, there is more than one way to linearize. The general idea behind linearizations, though, is to find a linear system that can be iterated such that, after many iterations, the solution to the linear system approaches the non-linear one.\n",
    "\n",
    "The system that we have just written down admits a particularly simple linearization. The non-linear system of equations we want to sovle (multiplying the equation above by $\\Delta x^2$ and doing some re-arranging) is\n",
    "\\begin{eqnarray}\n",
    "T_0 & = & 1 \\\\\n",
    "T_{i+1} - 2 T_i + T_{i-1} + \\frac{K'_i}{4 K_i} \\left(T_{i+1}-T_{i-1}\\right)^2 & = & 0 \\\\\n",
    "T_{N-1} & = & T_f\n",
    "\\end{eqnarray}\n",
    "Let us denote $T_i^{(m)}$ as the $m$th guess at the correct solution. Now consider the closely-related system of equations\n",
    "\\begin{eqnarray}\n",
    "T_0^{(m+1)} & = & 1 \\\\\n",
    "T_{i+1}^{(m+1)} - 2 T_i^{(m+1)} + T_{i-1}^{(m+1)} + \n",
    "\\frac{(K'_i)^{(m)}}{4 K_i^{(m)}} \n",
    "\\left(T_{i+1}^{(m)}-T_{i-1}^{(m)}\\right)\n",
    "\\left(T_{i+1}^{(m+1)}-T_{i-1}^{(m+1)}\\right)\n",
    "& = & 0 \\\\\n",
    "T_{N-1}^{(m+1)} & = & T_f\n",
    "\\end{eqnarray}\n",
    "The middle of the three equations can be rewritten simply as\n",
    "\\begin{equation}\n",
    "\\left(1 + f_i^{(m)}\\right) T_{i+1}^{(m+1)} - 2 T_i^{(m+1)} + \\left(1 - f_i^{(m)}\\right) T_{i-1}^{(m+1)} = 0\n",
    "\\end{equation}\n",
    "with\n",
    "\\begin{equation}\n",
    "f_i^{(m)} = \\frac{(K_i')^{(m)}}{4 K_i^{(m)}} \\left(T_{i+1}^{(m)} - T_{i-1}^{(m)}\\right).\n",
    "\\end{equation}\n",
    "It is clear that this system of equations is linear in $T_i^{(m+1)}$, and it can be solved for $T_i^{(m+1)}$ given any vector of input values $T_i^{(m)}$. However, it is also clear that if $T_i^{(m+1)} = T_i^{(m)}$ that it is also a valid solution to the original non-linear system. Thus we simply need to iterate our guesses for $T_i$ until they stop changing, i.e., until $T_i^{(m+1)} \\approx T_i^{(m)}$. This is a linearized system, and the process of iterating until $T_i$ stops changing is called Picard iteration.\n",
    "\n",
    "For the purposes of numerical solution, it is convenient to write the linear system of equations as a matrix equation. In fact, it is a particularly simple matrix equation:\n",
    "\\begin{equation}\n",
    "\\left[\n",
    "\\begin{array}{ccccc}\n",
    "1 & & & & \\\\\n",
    "1 - f_1^{(m)} & -2 & 1 + f_1^{(m)} & &  \\\\\n",
    "  & \\ddots & \\ddots & \\ddots & \\\\\n",
    "  &        & 1 - f_{N-2}^{(m)} & -2 & 1 + f_{N-2}^{(m)} \\\\\n",
    "  &        &        &        &  1\n",
    "\\end{array}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "T_0^{(m+1)} \\\\\n",
    "T_1^{(m+1)} \\\\\n",
    "\\vdots \\\\\n",
    "T_{N-2}^{(m+1)} \\\\\n",
    "T_{N-1}^{(m+1)}\n",
    "\\end{array}\n",
    "\\right] =\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "1 \\\\\n",
    "0 \\\\\n",
    "\\vdots \\\\\n",
    "0 \\\\\n",
    "T_f\n",
    "\\end{array}\n",
    "\\right]\n",
    "\\end{equation}\n",
    "Thus the problem is just a tridiagonal matrix equation (i.e., one where the only non-zero elements in the matrix are in the diagonal or just above or below it), which is easy to solve by standard methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Putting it Together\n",
    "\n",
    "We are now in a position to build a program to solve our example probem. We will work from the innermost loop outward. First, let us suppose that we have some initial guess for $T_i$, which we will call $T_i^{(0)}$. The first step is to solve for $T_i^{(1)}$, the next guess, by solving the band diagonal system. We then compare the two solutions, and, if they are close enough, stop iterating. If not, we repeat the procedure to generate new guess $T_i^{(2)}$, and so forth, until a pair of iterates are close enough. Then we're done. Now we can code this up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import the scipy library function we'll be using\n",
    "from scipy.linalg import solve_banded\n",
    "\n",
    "# Inner step: given a guess for T_i^(n), this function produces the next guess T_i^(n+1)\n",
    "# The arguments K and Kp are callables that return K and K' when given T\n",
    "def next_iter(T, K, Kp):\n",
    "    \n",
    "    # Construct f, as defined above\n",
    "    f = Kp(T[1:-1]) / (4.0*K(T[1:-1])) * (T[2:]-T[:-2])\n",
    "    \n",
    "    # Store the tridiagonal matrix in the compressed form used by solve_banded.\n",
    "    # The matrix is stored as a (3,N-2) matrix, with the upper diagonal in the\n",
    "    # first row (with a leading zero), the diagonal in the second row, and\n",
    "    # the lower diagonal in the third row (with a trailing zero)\n",
    "    A = np.zeros((3,T.size))\n",
    "    A[0,2:] = 1.0 + f   # Upper diagonal\n",
    "    A[1,0] = 1.0        # Upper left corner of diagonal\n",
    "    A[1,1:-1] = -2.0    # Diagonal\n",
    "    A[1,-1] = 1.0       # Lower right corner of diagonal\n",
    "    A[2,:-2] = 1.0 - f  # Lower diagonal\n",
    "    \n",
    "    # Right hand side matrix\n",
    "    b = np.zeros(T.size)\n",
    "    b[0] = T[0]\n",
    "    b[-1] = T[-1]\n",
    "    \n",
    "    # Solve the tridiagonal matrix equation using solve_banded\n",
    "    T_new = solve_banded((1,1), A, b)\n",
    "    return T_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Outer step: this function takes a value of Tf, and iterates until T converges\n",
    "# The optional argument T0 is an initial guess\n",
    "def Tsolve(Tf, N, K, Kp, tol=1.0e-6, T0=None):\n",
    "\n",
    "    # If not given an initial guess, use a linear ramp; this is the exact solution\n",
    "    # for constant K\n",
    "    if T0 is None:\n",
    "        T = np.linspace(1, Tf, N)\n",
    "    else:\n",
    "        T = T0\n",
    "        T[0] = 1.0\n",
    "        T[-1] = Tf\n",
    "    \n",
    "    # Iterate\n",
    "    itcount = 0\n",
    "    while True:\n",
    "        T_last = T    # Save last state\n",
    "        T = next_iter(T, K, Kp) # Iterate\n",
    "        itcount = itcount + 1\n",
    "        # Check for convergence\n",
    "        if np.amax(np.abs(T - T_last)) < tol:\n",
    "            break\n",
    "            \n",
    "    # Return converged solution and iteration count\n",
    "    return T, itcount"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define an example K and K' function pair: K(T) = T^2, K'(T) = 2 T\n",
    "def K(T):\n",
    "    return T**2\n",
    "def Kp(T):\n",
    "    return 2.0*T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example application: solve for T at a range of Tf values; keep track\n",
    "# of iteration count\n",
    "Tf = np.arange(2, 10.0001, 0.1)\n",
    "N = 1000\n",
    "Tsol = np.zeros((Tf.size,N))\n",
    "it = np.zeros(Tf.size, dtype=int)\n",
    "for i in range(Tf.size):\n",
    "    T, itcount = Tsolve(Tf[i], N, K, Kp)\n",
    "    Tsol[i,:] = T\n",
    "    it[i] = itcount"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 10)"
      ]
     },
     "execution_count": 120,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbwAAAFBCAYAAAAMvBoYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl02+d95/v3DytBbCQIriBBgKskitRG7ZJtWZbXtNnb\nZmuSpk2bNl3vtJMznbmZtve0vWe63M7c3k5zkzRt2kyTprnN5ji1E1te5EimFooSKe4kuBMgQJDY\ngd/vuX+Aou0kjiWLJLg8r3N0aFkg8ZVN4oPvsypCCCRJkiRpu9MVugBJkiRJ2ggy8CRJkqQdQQae\nJEmStCPIwJMkSZJ2BBl4kiRJ0o4gA0+SJEnaEdY98BRF+byiKPOKotx4zb9zKYrytKIogysfS9e7\nDkmSJGln24gO7wvAoz/07z4FfE8I0Qx8b+X3kiRJkrRulI3YeK4oig/4lhBi78rv+4EHhBAziqJU\nA88JIVrXvRBJkiRpxyrUHF6lEGIGYOVjRYHqkCRJknYIQ6ELeDOKonwc+DiA1Wo9tGvXrgJXJEmS\nJN0pVUui5uLk1BhqLgGoAOh0RZjNlRgMjnv6+pcvXw4JIcrv5LGFCrw5RVGqXzOkOf9GDxRCfAb4\nDEBnZ6fo6uraqBolSZKkuyCEIJkcIxx5mXD4JSKRH5DLpQA9xcUHKS09RmnpcUpLjmAyla3JcyqK\nMn6njy1U4H0D+DDwpysfv16gOiRJkqR7kE4HiUQuEI5cIBx+iXR6BgCzuYpy91lKXSdwlR7HbK4s\ncKUbEHiKovwv4AHArSjKJPBp8kH3FUVRPgYEgPeudx2SJEnSvcvllolELhKOXCASuUA8PgiAweCk\ntPQYrtJP4HKdwGLxoShKgat9vXUPPCHE+97gj86u93NLkiRJ90bT0kSjVwmHXyIceZnl5esIoaLT\nFVHi7KS66p2Ulp7Abt+DougLXe5PtOkXrUiSJEkbRwiV5eXefAcXvsBitAtNS6Eoeuz2Duq9v4zL\ndRKn8wA6nbnQ5d4VGXiSJEk7XDI5RTj8AuHwi4QjF8jlogBYrc3U1PwsLtdJSkuOYDDYC1zpvZGB\nJ0mStMPkcnEiiz9YDblEYhS4vdDkoXzAlR7HbN5eW6Rl4EmSJG1zQmgsL98gHH6RhfCLRKNXECKL\nTldEaelRPJ4PUOY6TXFx46ZbaLKWZOBJkiRtQ6nUTH6hSfgFwpGXyGYjANhse/DW/QIu1ylKSg5t\nuXm4eyEDT5IkaRtQ1SSLi5dYCL9IOPzC6nYBk6mcsrIHKHOdptR1ErPJXeBKC0cGniRJ0hYkhCAW\n7ye88DwL4RdYXOxCiAw6nYkS5xGqq9+Ny3Uam7V1Ww9T3g0ZeJIkSVtELrdMOPwSCwvnWQg/Tzo9\nC4DV2kJd7YdwuU5TUnIYvb6owJVuTjLwJEmSNqnbXdzCwnkWFp5bWWySw2Cw4yo9RVnZ/bjKTlNk\nrip0qVuCDDxJkqRN5I26OJttN17vL1JW9gBOx350OmOBK916ZOBJkiQV0Ou7uPNEo5d/pIsrK7tv\nUxy+vNXJwJMkSdpgsosrDBl4kiRJGyCRGCUU+j6h0PdZjHbJLq4AZOBJkiStA03LshjtYiH0LKGF\n768e32WztsourkBk4EmSJK2RbDZCaOE8odD3CYefJ5dbRlFMuEqPUVv7YdxlZ7BYagtd5o4lA0+S\nJOktEkIQTwytDlVGo1cADZPJTUX5Y7jdZygtPYnBYC10qRIy8CRJku6KpqWJLL5CKPQ9QqFnSaUm\nALDb2vD7fg23+0Hs9r0oiq7AlUo/TAaeJEnSm8hmI4RC3ycY+h7h8IuoahydrgiX6yS++l+mzH1G\nbv7eAmTgSZIk/RjJ5ATB0DMEg0+zuPgKoGE2V1FV+dO43WcpLT2GXm8pdJnSXZCBJ0mSxMoG8Fgv\nweDTBEPPEIv1AflVlT7fJyh3n1sZqpQHMW9VMvAkSdqxNC3L4uIrBEPPEAo+TSo9DegocR6iuek/\n4XY/RHFxfaHLlNaIDDxJknaUXC5OOPwCwdDThELPkstF0enMuFyn8ft/E7f7DCZTWaHL3LZisRiB\nQIBAIEBbWxt1dXUb9twy8CRJ2vYymVB+0UnwacKRF9G0DAZDCeXus5SXn8PlOoVeX1zoMrcdIQQL\nCwurARcIBAiHwwAYDAbcbrcMPEmSpHuVSs8SDP478/NPrS46KSqqxVPzfsrLz+F0dqLTyZfAtZTL\n5ZiZmSEQCDAxMUEgECCRSABgsVjwer0cOnQIr9dLdXU1BsPG/veX/7clSdo2kslJgsHvMj//HaJL\nVwGwWpvx+X6VivJHsNl2y0UnayiVSq0GWyAQYGpqilwuB0BpaSnNzc14vV68Xi9ut7vg/+1l4EmS\ntKUlEqPMzz/FfPAplpdvAGCz7aGh4XeoKH8Uq7WxwBVuH9Fo9HXDk3NzcwAoikJ1dTWdnZ14vV7q\n6uqw2+0FrvZHycCTJGlLEUIQjw8wH/wuwfmniMX7AXA49tPU+B+pqHgUi8Vb4Cq3vtvzb+Pj46u/\notEoACaTidraWh544AG8Xi8ejwez2Vzgit+cDDxJkjY9IQTLsZsEVzq5/M0DCiXOTpqb/zMV5Y9Q\nVFRT6DK3NE3TCAaDrwu4WCwGgNVqxev1cvz4cbxeL5WVlej1+gJXfPdk4EmStCnd3gg+N/dt5uaf\nJJWaQFH0lJQcpa72o5SXP4zZXF7oMrcsVVWZm5tjbGyM8fFxAoEAyWQSAIfDgd/vp76+nvr6+k0x\n/7YWZOBJkrSpxGL9zM1/m7m5b5NMjqEoBlylJ1YOZj6LyeQqdIlbUi6XY3p6erV7CwQCZDIZIL/A\npLW1FZ/PR319PSUlJdsi4H6YDDxJkgouHh9hbv7bzM9/m3h8ENBRWnqMeu8vUVHxCEZjaaFL3HKy\n2SyTk5OrATcxMbG6grK8vJyOjo7VDs7hcBS42o0hA0+SpIJIJgMrw5XfXjm3UqGk5DCtLX9AecWj\nmE3uQpe4pWSzWSYmJhgdHWVsbIypqSk0TQOgqqqKQ4cOrQac1boz7+eTgSdJ0oZJpaZXhyuXl3sA\ncDgO5BeeVDwmr9i5C7lcjsnJScbGxhgdHWVychJVVVEUhZqaGo4dO4bP56Ourg6LRd7qADLwJEla\nZ5nMAnPzTzI3942VG8HBbm+nqelTVJQ/jsXiKXCFW4OqqkxPT692cIFAYHWIsrq6miNHjuD3+/F6\nvRQVFRW42s1JBp4kSWtOVRMEg08zO/d1wuEXEULFZm2lseE/UFHxuLyB4A5omsbMzMzrAu72IpOK\nigoOHTqEz+fD5/PJDu4OycCTJGlNaFqWcPhFZue+QTD4NJqWpMhcg9f7S1RV/jQ2W2uhS9zUNE1b\n3SYwOjrK+Pg46XQaALfbTUdHB36/H5/Pt2Pn4O6VDDxJkt4yIQRLS1eZnf0Gc/PfJpsNYzA4qa56\nB5VVb6fEeQhF0RW6zE3p9kkmIyMjq13c7X1wLpeLtra21YDbjMd0bUUy8CRJumvx+DCzc19nbvab\nJFMBdDozbvdZqirfTlnZfeh0pkKXuCnFYjFGRkZWfy0tLQHgdDpX98H5/X6cTmeBK92eZOBJknRH\n0pkQc7PfYHbu31hevgno8hvC/Z+kvPxhDAbZhfywdDpNIBBgeHiYkZER5ufnASgqKsLv93P69Gka\nGxspLS3dlhu9NxsZeJIkvSFNSxMKPcvMzL+yED6PECp2ezvNzf+ZyoonMJsrCl3ipnJ7JeXtDm5i\nYgJN09Dr9Xi9Xs6ePUtDQwPV1dXodHKod6PJwJMk6XWEECwtX2dm5mvMzX2TXC6K2VSJt+4Xqa5+\nF1ZrU6FL3DSEEIRCodWAGxsbW11oUl1dzfHjx2loaMDr9WI0GgtcrSQDT5IkIH9D+Ozs15mZ+RqJ\nxBA6nZny8keornoXLtcJFGXrnY6/HmKx2OoQ5cjICMvLy0D+PMq9e/fS0NAgV1JuUjLwJGkHU9Uk\nweDTzMx+jXD4JUDD6exk164/prLicTkvR36YcmJigqGhIYaHh5mZmQHAYrHg9/tpbGzE7/fjcslD\nrTc7GXiStMMIIYhGLzMz86/MzT+JqsYoKvLg8/0q1VXvpLjYV+gSCy4SiTA0NMTQ0BCjo6NkMhkU\nRaGuro4zZ87Q1NQk5+G2IBl4krRDZDIhZmb/jenpr5BIDKPXF1NR/hjV1e+ipOTIjt4vl8lkGBsb\nWw25cDgM5LcLtLe309TUhN/vl0d2bXEy8CRpGxNCJRx+kanprxAKPYMQOZzOg+ze9X9SUfEYBsPO\nnGcSQjA3N8fw8DBDQ0MEAgFUVcVgMOD3+zly5AhNTU2UlZXJ7QLbSEEDT1GU3wZ+ERBAD/BRIUSq\nkDVJ0naQTE4yM/NVpme+Sjo9g9Hooq72w9TU/MyOXWWZTCZXA25oaIhYLAbkz6U8evQojY2NcjXl\nNlewwFMUxQP8BrBHCJFUFOUrwM8BXyhUTZK0lWlammDwGaanv0I48hIALtcpmpt/n3L32R13+okQ\ngvn5eQYGBhgcHGRiYgIhBBaLhYaGBpqammhsbNwxl59KhR/SNAAWRVGyQDEwXeB6JGnLicUGmJ75\nCrOz/0Y2G6HIXIPf/xvUVL+HoqKaQpe3oTKZDKOjo6shd/vorqqqKk6fPk1zczMej0cuNikgoQly\nwQSZyRimOjvGiuINe+6CBZ4QYkpRlD8DAkAS+HchxL//8OMURfk48HEAr9e7sUVK0ialaWnm57/L\n1NSXWIy+gqIYKS8/R031e3G5Tu6oPXPhcJjBwUEGBgYYGxtDVVVMJhMNDQ088MADNDU1yS6uQIQQ\nqAspMpPLZCZjZCaXyU7HEJn8TezOx/wbGniKEGLDnux1T6wopcC/Aj8LLAL/AnxVCPGPb/Q5nZ2d\noqura4MqlKTNJ5kMMDX1v5ie+SrZbBiLxYun5n1UV78bk6ms0OVtiFwuRyAQWA25hYUFAMrKymhu\nbqalpQWv14vBUOgBrJ1HXc7kw20iH3DZyWW0RP6SWgw6TDVWTLV2jLU2TLV2DG4Liu7eFgUpinJZ\nCNF5J48t5HfEQ8CoECIIoCjK14ATwBsGniTtRJqWY2HhWSan/olw+AUURY+77EE8ng+sdHPbf3gu\nFoutDlMODw+TyWTQ6/X4fD4OHz5Mc3MzZWU7I/A3Cy2tkp1aJjMRWw05dTF/rBoKGCutWNrcGOvy\n4WasLEbRF/Z7tZCBFwCOKYpSTH5I8ywg2zdJWpFOzzE1/RWmp/+ZdHoWs6kSv+83qKn5GYqKqgtd\n3roSQhAMBunv76e/v5/JyUkAHA4H7e3tNDc309DQgMm0sxbiFIpQNbKziVe7t4llcvOJ/Pp6QO8q\nwuS1YzpZk5+Xq7GhM22+YfVCzuFdVBTlq8AVIAdcBT5TqHokaTMQQiMcucDU1JdW9s2puFynaWn5\n33GXnUWn277DdKqqEggEVkMuEokAUFNTw5kzZ2htbaWyslLui1tnr5t3Wwm3zHQccvl5N53VgKnW\nTnG7G2OtHVOtDb1ta7zxKOhPjxDi08CnC1mDJG0GuVyMmdn/j8nJfyCRGMFoLKWu7hfw1Pzctj7q\nK5VKMTw8TH9/PwMDA6RSKfR6PQ0NDZw8eZKWlha54GSdqfFsPtQCS/mFJRPLiGR+3k0x6jB6bNiO\nVWOqs2Oqs6MvNW/ZNx3b9+2iJG0BicQYk5NfZHrmq6hqDIe9gz27/4yKisfR682FLm9dLC4uMjAw\nQH9/P6Ojo2iahsViobW1ldbWVhobGzGbt+ffvdCEKsjOxvPhFlgZmgwl83+oy8+7Fbe783NuK1sG\nFP3WDLcfRwaeJG0wITTC4ReZmPx7FhaeQ1EMVFQ8Tl3th3E69xe6vDV3+xivvr4++vv7mZ2dBfKr\nKo8dO0Zrayt1dXVyb9w6UJczq+GWDiyRnYwhsitDkzYjJq8D6+FKTHUOjLWbc95tLcnAk6QNkh+2\n/NrKsOUoJpMbv+838Hjet+1uDtc0jampKfr6+ujr6yMSiazeNnDu3DlaW1txu92FLnNbETmN7Eyc\n9O3uLbCEGllZNalXMNbYsB6uwlRvx1Tn2NJDk2+VDDxJWmeJxCgTk19kZuZf88OWjn207fkLKioe\nRafbPkN3qqoyPj6+GnKxWAydTkdDQwOnTp2itbUVm81W6DK3jVw0ne/exlcWlkwtQy6/bFLvNGHy\nOjCdqMl/rLGhGGUHLQNPktaBEILFxVcITHyOUOh7KIqByorHqa39+W01bJnNZhkeHl4drkylUhiN\nRpqamti9ezctLS3ySp01IFSN7HSc9Fg0372NL6EuZfJ/aFAweezYjtfktwZ4HRic2+eN1FqSgSdJ\na0jTcswHv0Mg8DmWl3swGkvx+X6NWs8HMZvLC13emkilUgwODtLX18fg4CDZbJaioiJaW1vZtWsX\njY2Ncn/cPdISWdIrwZYeWyI7ubw696YvMWPyOzF57Zi9DozVVhSD7N7uhAw8SVoDudwy09NfYWLi\nC6TS01gsPlpb/4jqqnei11sKXd49SyaT3Lp1i97eXkZGRlBVFavVyr59+9i9ezc+nw+9fnsveFgv\nQgjUcIr02FI+4MaXyM0l8n+oIz/3dqQKU70Ds8+B3iG7t7dKBp4k3YNUapqJiS8wNf1lVDVGSckR\nWlo+jdv94JY/8iuZTNLf38/NmzcZHh5G0zScTidHjhxh9+7d1NbWypWVb8Grw5NLZMajpMeX0Jaz\nAChmPaZ6B8Ud5Zh8Dkx19m2/cnIjycCTpLdgafkGgcBnmZ9/EoCK8sfwej+Gw9FR4MruTSqVel3I\nqaqK0+nk2LFj7NmzB4/Hs+NW9t2rnzg8WWqmqLEkH271zvx5k/d4mLL0xmTgSdIdEkIQibzM+Pj/\nJBx5Cb3eRl3tR6it/TAWi6fQ5b1lqVSKgYEBbt68ydDQEKqq4nA4OHLkCG1tbTLk7pK6lCY9ukR6\nNEp6NPrqmZNyeLLgZOBJ0psQQiMYeprx8b9laakbk6mcpsbfw+N5PwaDvdDlvSXpdHo15AYHB1FV\nFbvdzuHDh1dDTg5XvjkhBGokvRpumdEouYUUAIpJj6nenh+erF8ZnjTL4clCkoEnSW9A07LMzX2D\nsfHPkEgMYSnyrixEefeWPPYrm80yODhIT08Pg4OD5HI57HY7nZ2dtLW1yTm5OyCEIDefWO3gMmNR\n1Gh+e4Cu2IDJ58R6rBqz34mx2ratjuXaDmTgSdIPUdUk09NfZjzwWdLpGWy2XbTt+UsqKh7fcrcV\naJrG6OgoPT099PX1kU6nsVqtHDx4kLa2Nnmk15sQmsifXnK7gxuLosXzByvr7CbMfgdmvxOz34mh\nQs6/bXZb66dXktZRNrvE5OQ/MDH592SzYZzOTna1/hFlZQ9sqTksIQRTU1P09PRw48YN4vE4ZrOZ\n3bt3097eLrcQ/AQip5GZiq0OT6bHlhBpFcjf+VbU6loNOH1Z0Zb6vpBk4EkS2WyUiYm/Y2LyC+Ry\ny5SVncFX/yuUlHQWurS7EgwG6enpoaenh0gkgl6vp6WlZfXCVKPRWOgSNx2hamQmY6RHoqRHFsmM\nLa2uoDRUWCjeX47Z78Tkc2Io2XrD2NLrycCTdqxsNkIg8HkmJv8BVY1RXv4Ift8nsdv3FLq0OxaN\nRrlx4wY9PT3Mzs6iKAp+v5/77ruPXbt2YbFs/U3va0loguxUjPTIIqnhaD7gMvkOzlBZnD9c2e/E\n7HdsmUtNpTsnA0/acTKZMIGJzzM5+Q+oapyK8sfw+T+J3bar0KXdkXQ6TV9fH93d3YyOjgLg8Xh4\n9NFHaWtrw27fmitH18PqHNzIIunh/Dzc7SFKQ7mF4oMVmBucmBucMuB2ABl40o6RySwQCHyOyakv\noqpJKioew+/7JDZba6FLe1OapjE2NkZ3dze9vb1ks1lKS0u5//776ejooKysrNAlbgpCy6+iTA2/\nJuBWbu82lBVRvK98JeBK0DtkwO00MvCkbS+TWWA88P8yNfVPqGqSysq34fP9GjZrc6FLe1OhUIju\n7m66u7tZWlrCbDbT3t7Ovn378Hq9O37RhBCCXChJemhxdR7u9ipKvasIy54yzI1OzI0l8gYBSQae\ntH3lcsuMBz7LxMTfoapJqip/Cp/v17BaGwtd2k+USCS4ceMG3d3dTE1NoSgKjY2NnDt3jl27du34\nxSfqcob08CKpwUXSQ4uo0fwlp3qnOb+KsqEEc4MTg0teSyS9ngw8adtR1QQTk19kfPxvyeWiVFQ8\nToP/tzZ10KmqytDQENeuXWNgYABVVamoqODhhx+mvb19R8/LaRk1vw9uJeCys3EAFIuBokYn5jN1\nmJtKMMhtAltGLpNhYTKAtdSFrdS1Yc8rA0/aNjQtzdT0lxkb+2symRBlZQ/Q2PA72O1thS7tDS0s\nLHDlyhW6u7uJxWIUFxdz+PBh9u3bR1VV1Y58AReqIDO1THpwkdTQIpnAEqgC9ApmnwPHIz6Kmkow\nemxyo/cmJ4QgFlkgND7G/PgoocAYwfFRwtOTCE3jzEd+mYOP/dSG1SMDT9ryNC3H7Oy/MTr6V6TS\n05SUHKV9719v2n10mUyG3t5erl69yvj4OIqi0NLSwoEDB2hubt5xm8JfOw+XGlwkPbKISOVXUhpr\nrNhOeihqyt8oIK/K2bxud23B8VGCK8EWDIyRWl5afYzdXU55vZ+mw8cpr/dT07qxK6Nl4ElblhCC\n+eBTjIz8BYnECA57B7t2/wmu0pObrjMSQjA9Pc3Vq1fp6ekhnU7jcrk4e/Ys+/fv33FDlloiS2po\nkdRAhPTga+bhSswUt5djbirB3Ci3CmxGt7u24PgowfF8sIUCY6tdG4DBZMbtraf58DHK6/2Ue/24\n630UWW0FrV0GnrQlLS52MTj0pywtXcVqbaaj/W9wu89tuqBLJBL09PRw5coV5ubmMBgM7Nmzh4MH\nD1JfX7/p6l0vQhNkJpdJD0RIDUTITCyDAKVIT1FjCeYzdRQ1l6B3yXm4zUTNZVmYnGB+bITg2Mhq\n55aKLa8+xlFegdvro/nIcdxeP+X1fkqqqtDpNl83LgNP2lLi8RGGR/4bweC/YzZVsnvXn1Jd/S4U\nZfP8cAkhGB8fp6uri76+PlRVpbq6mieeeIK9e/fumNNP1KUMqYEIqcEI6cEIWiIHChg9Nuxn6ihq\ndWGqtcsbBTaJdCJBcHyE+bFR5seGmR8bYWEigKau7GO83bUdPUG510d5vR+3t/Bd292QgSdtCelM\niNHR/8709D+j01loaPgdvHW/gF6/ecIjmUzS3d1NV1cXoVAIs9nMoUOHOHDgANXV1YUub92JnEZ6\nfCk/TDkQITuTX02psxkp2uWiqKUUc3MpeuvO3lZRaEII4osR5seGCY6NMj+aD7fFuZnVx1gcTip8\nDfieeDsVvgYq/I2UVFVvyq7tbsjAkzY1VU0QCHye8cBn0LQ0npr34/d/EpPJXejSVk1NTdHV1UVP\nTw+5XA6Px8Pb3/522traMJm29xxUbiGZ7+IGIqSHFxEZLb+ast6B41EfRS2lGKuscjVlgQhNIzI7\nsxJuI8yv/EpEF1cf46ysosLXQNv9Z6nwN1Lha8Ba6tqWQ8sy8KRNSQiNmdmvMTL8F6Qzc5SXP0JT\n4+9SXOwvdGlAfqXljRs3eOWVV5iZmcFoNNLR0UFnZyc1NTWFLm/dCFUjPbZE6laY1K0wuWASyJ9q\nUnywMt/FNTrRmeVLy0a7vUpybnQ4PzQ5OkJwfJRsOn8Du06vp6zWi39/JxX+BirqGyj3+TEXWwtc\n+caR35XSphONXqF/4A9ZXu7B4djP3r3/fdNsMZifn6erq4vu7m7S6TTl5eU8/vjjdHR0UFS0PU/2\nUGMZUv2RfMgNRPKHL+sVzA1OrEerKdrlkpu+N1g2kyY0PsbcyBBzo0PMjQyxMBlAU/PbOUwWC+X1\nfvaeOUeFr4FyXwNltV4MO/yUHhl40qaRSs8yPPTfmJ37N8ymStr2/AWVlT+FohT2Rm5N0+jv7+fi\nxYuMjY2h1+vZs2cPnZ2d2/I8SyEE2en4aheXmcyvqNTZTVja3Vh2uzA3laIzb+35nK0im0kTHBtd\nDbb5kSFCk4HVLQAWu4PKhib8Bzqp9DdS7mugpKIKRd5k/yNk4EkFp6opAhOfY2zsbwAVX/2vUl//\nKxgMhR1qSSQSXL16lUuXLhGNRnE6nTz00EMcOHAAq3V7DQNpaZX0UITUrQjJ/jDaUia/orLWjuOh\neop2uTBWy7m49ZZNpwiOj+Y7t5Fh5kbzndtquDmcVDY00XDoKJUNjVQ2NGEvK992b7rWiww8qWCE\nEASD32Vw6E9IpSYpL3+U5qZPYbHUFbSuubk5Ll26RHd3N7lcDp/Px6OPPkpLS8u2OgUlt5gm1bdA\nsneB9EgUVIFi1lPUUppfVdlaKjd+r6MfCbeRQRamJn4k3Jo6j1LR0ESlvwl7mVuG2z2QgScVRDw+\nTP/Ap4lEXsZmbeXAgX/EVXq8YPVomsbAwAAXL15kdHQUg8FAR0cHR44coaqqqmB1rSUhBNnZBKmb\nIZJ9YbJTMQAMbgu24zUU7XZhrnegGORQ2FrLZTIEx0eZHR5YCbghFiYnECIfbsXOEir9jTQdPibD\nbR3JwJM2lKomGRv7a8YDn0Wvt9Da8gfU1PwcOl1hvhVTqRRXrlzh0qVLLC4u4nA4eOihhzh48CDF\nxcUFqWktCVWQHouS6s13cmokDQqY6uw4HvVh2VOGsWLr/z03E01TiUxPMTM0wOzQALPDgwTHR1c3\ncBc7S/Kd25HjVPqbqGxowuYqk+G2AWTgSRsmFPo+/QN/QCo1SVXVO2lq+hTmAu2nW1xc5OLFi1y+\nfJlMJkN9fT0PP/wwra2tW37YUkur+b1xvQskb4XzN34bFIqaSnGc8VK024XeLocq14IQguWFILPD\ng/lwGxpgdmSIbCq/XcNksVDV2Myht72D6sYWKhubZedWQDLwpHWXTE4xMPiHhELPYLU2c/DAlygt\nPVqQWqanp7lw4QI3b94EYO/evRw/fnzL751TlzMkexdI9S6QGl6EnEBXbMCy25W/9btZrqpcC8nY\nMnNDA8zs7mAlAAAgAElEQVQMD6yG3O1N3HqDgXJfA233P0hVYwtVjS24ajxyteQmIgNPWjealiEQ\n+DyjY/8DUGhq/I/U1X0UnW5j9wJpmsbQ0BAXLlxgbGwMk8nEsWPHOHr0KCUlJRtay1rKLaZJ3giR\nvBEiM74EIr8B3Ha0mqI9ZZh9TnlO5T3IplPMj47kg204372tHr+lKLhqavHvP7QSbs246/07fp/b\nZicDT1oX0aVu+vo+RTw+QHn5w7Q0/xeKija2i8pms/T09HDhwgVCoRAOh4OHH36YgwcPbtlN4rmF\nJMkbC/mQm8ifWG+sKsZx1otlrxtDZbEcLnsLhBBE52aZGbzF9OAtZgb7V+bd8hu57WXlVDU10372\nEaoaW6hsaMK8DeZ4dxoZeNKaUtUEwyN/ycTEFzCbK+jo+Azl7rMbWkMqlaKrq4uXX36ZeDxOVVUV\n73rXu2hra9uS83PZ+cRqJ5edzh/IbPTYcDziw7K3DGO5fOG9W5lUktmhwdcFXHIpCoCxyEJ1UzOH\nf/rdVDW1Ut3UgrWktMAVS2tBBp60ZsLhl+i79fukUhN4PB+gqfF3MRg27mLTeDzOxYsXuXTpEqlU\nioaGBk6dOoXf799SXc/t7QPJGyGSPSFy8wkATF47zif8WNrcGFxbs0MtBCEEkZlpZgZv5QNu4Bah\nwPjqloDSmloaDhymurmVmpZdlNV5t/ytANKPJwNPumfZ7CKDg3/MzOy/Ulzs5+DBf6a05PCGPf/S\n0hIXLlzg8uXLZLNZdu/ezalTp/B4PBtWw1rIzsZJdAfzIRdKggJmvxPbsUYsbWXoneZCl7glpBMJ\nZocGmB7sY2awn5nB/tULS02WYqqbWzn6rmPUNLdS1dyKxbazbpvfyWTgSfdkbv47DAz8V7LZCL76\nT+Dz/Tp6/ca8MC8sLPDSSy9x7do1hBB0dHRw8uRJKioqNuT510I2mCB5PUSiO5jv5BQwN5ZgO+3J\nh5w86eQnuj33NtXfy1R/LzMDtwhNBkAIUBTKPHU0HT5OTcsuqptbKfPUyVWTO5gMPOktyWYX6e//\nNHPz38Jub2P/vr/Dbt+zIc89NzfHCy+8wM2bN9HpdBw6dIgTJ05QWro15llykRTJ60ES10P5004U\nMPkclLyjEctetwy5n0DN5QiOjawG3HR/H/HFCABmq5Xq5l20HDtFdcsuqptadtTVN9Kbk4En3bVQ\n6Fn6bv0nstkwDf7fpr7+VzbkpJS5uTnOnz9Pb28vJpOJEydOcOzYMez2zT8kpS5lSPQESXYHyQRW\nhtfq7DifaMDS4cYghyt/rHQizvTALab7e5nq72NmqJ9cOg2As6ISb/t+PK178LTupqzWK7s36Scq\naOApilICfBbYCwjgF4QQLxeyJumN5XIxBof+mOnpL2O1trB/3+c2pKubn5/n/Pnz3Lx5E5PJxH33\n3cexY8c2/dFfWiJLoidEsjtIejQKAozVVhyP+ihud2MosxS6xE1nKTTP1K18uE3fuklwYhyEQFF0\nVPgbaH/wYTytbXhad2NzlRW6XGmLKXSH91fAU0KI9yiKYgI29yvYDhaJXKS37/dIpaap9/4yDQ2/\niU63vl3JDwfd6dOnOX78+KYOOpHTSN0KE78yT6o/DKrAUG7J75PrKJfnVr6G0DRCE+NM9t1gqr+P\nqf5eYgshIL81oKZlF8ePnMCzaw/Vza2YiuQbBOneFCzwFEVxAPcBHwEQQmSATKHqkX48TUszPPzn\nBCY+j8Xi5dChf6bEeWhdnzMYDHL+/Hlu3LiB0Wjk1KlTnDhxYtMGndAEmfElElfnSVwPIVI5dDYj\ntuM1FB+owFhj3VLbItaLpqoEx0eZ6O3Jh1zfTVLx/I0NNlcZntY91LTuwbNrD+VeH7otuGdSunNq\nVkMTAqNp4/4/F7LDawCCwN8pirIPuAz8phAiXsCapNeIx0e4efO3WI7dxOP5AM1Nn0KvX7/QiUQi\nPPvss1y/fh2j0cjJkyc5ceLEpr1sNTufyIfctXnUSBrFqMOy103xgQrMjSU7/lgvNZdjbmSQyb6b\nTPb2MNXfSyaZP1S5pLKapiPHqd29l9rde3GUV8g3BdtUMpYhMptgcTZBZDZOZC5BZDbBcijJ/e9v\npe30xm0fKmTgGYCDwK8LIS4qivJXwKeA//LaBymK8nHg4wBer3fDi9yJhBDMzPwL/QN/iF5fREf7\n31Je/tC6PV8sFuP555+nq6sLnU7HiRMnOHny5KYMOjWWIdEdJHF1nuxkfoWluakEx8P5q3Z28gHN\nuUyG2aEBJvp6mOy7yfRA3+oCE5enjt2nHsCzey+1u9uwuwpzS4a0PjRNsLyQJDKbWAm3+Oo/p+LZ\n1cfpjTpKKoqp8NppOVJJRb1jQ+tUhBAb+oSrT6woVcAPhBC+ld+fBj4lhHjijT6ns7NTdHV1bVCF\nO1M2u8St/t9nfv5JSkuPs2fPn1FkXp8LUFOpFBcuXODll18ml8tx8OBB7r//fhyOjf0heDNCXZmX\n65rLz8tpYKyxUnygkuJ95egdO3MbQTadYnrgFpN9N5jsvcHMUD9qNguKQnldPbV72qnds5faXW0U\nO7fuId3SqzKpHIsrHVr+Yz7YFucTaLlXs8RiN1JaZaWkqpjSymJKq6yUVhVjcxWh061tJ68oymUh\nROedPLZgHZ4QYlZRlAlFUVqFEP3AWaC3UPVIsLjYxc2bv006M09j4+9R7/0lFGXtl3lns1m6urp4\n/vnnSSaTtLW1cebMGdzuzfWuPzsXJ/7KHIlr82ixLDq7EdvpWqwHKzBWbr7uc72puSwzg/0Eblxn\n4uZ1pgduoam5lRWUjex/5G3U7t6LZ9ceeXrJFiaEIL6YITIXXxmGzAfb4lyCWCS9+jhFp+Ast1BS\nWUz93jJKq/LBVlJZTJF1c94aUehVmr8O/NPKCs0R4KMFrmdHEkJjfPwzDI/8ORZLLYcOfQWnY9+a\nP4+maXR3d/Pss8+ytLREY2MjZ8+e3VR30WnJHInuIPHLc2QnlkGnYNntorizkqIW146al9M0lfmR\nYQI3rxO40c1Uf29+iFJRqPQ3ceiJt1O3p52a1j3y5oAtSGiCpYUUkZk44ZVfkZl8x5ZNq6uPMxbp\nKa2y4mktzYdaZb5zc5Zb0Bu21r7HggaeEOIacEetqLQ+stlFent/l9DC96moeJzdu/54XQ58HhkZ\n4bvf/S5zc3PU1NTwjne8g4aGhjV/nrdCaIL0yCLxrjmSNxYgp2GsKsb5RAPFB8p3zMknQtMITQaY\nuNFN4OZ1JntvkE7k15C56+ppf/BhvG37qN29lyKbrcDVSndK0wRLwWQ+0GZXwm06373lstrq46wl\nZkqritl1ojo/DFmdH4Ysdpi2zYKiQnd4UgEtLV2n58YnSafnaWn5NLWeD635N3YwGOTpp59mYGCA\nkpIS3vOe99DW1rYpfoDUaJr4K7PEu+ZQF9MoRQasnZVYOysxemybosb1JIRgcXaawI3rBG7mhylv\nX5FTUlVN6/HT1O3toG5Pu7weZwvQVI3o7WCbiROeSRCeyQebmns12GylZlzV+Y7NVW3FtRJs5uLN\nOQy5lmTg7UBCCKamvsTA4P+B2eTm0MF/xuncv6bPEY/HOX/+PK+88gomk4lz585x5MgRjAW+EVpo\ngtRghPjFWVK3FkDLr7J0PurD0laGYtzeqywTS1ECN7oZv36V8evXWF4IAvl9cP79h/Du3UddWzsO\n99Y5gHunUXMa0fnk64YhwzP5OTZNfXXhiL2sCFe1lbrdrtcFm8myc1/2d+7ffIdS1QR9t36fublv\nUFZ2P217/hyjce3evedyOS5dusT58+fJZDJ0dnbywAMPFHyLgbqcId41S/zSLGokjc5qxH66FuuR\nqm19xFcum2W6v4/xnquMX7/K3OgwCIHZasW7dx9H3/levHv3UVJVs+072q1G0wTR+QTh6TgL03HC\n0zHC03Gi80k0bSXYFHC4LbiqrfjayyhdDTYrxh28ReaNyMDbQZLJKa73/AqxWB8N/t/G5/vVNV2F\nOTAwwHe+8x0ikQjNzc2cO3euoFf1CE2QHl4kfmmW5M0F0ATmBifOR/35bm6LTbjfCSEE4akJxq9f\nZez6VSZ6e8il0+j0eqqbd3HyvR+gvuMAlY1N8pLTTSK/KjLNwnSchal8qC1MxYjMJlBvz7Ep4HRb\ncNVYadhf/ppgK8awgSeVbHUy8HaISOQiPTc+iRBZ9u37LO6yB9bsa4fDYZ566ikGBgZwu9188IMf\npKmpac2+/t1S41kSXbPELs2iLqTQFRuwnazBeqQKY/n2W02YWIoy3nNtZZjyKrHwAgCl1R72PnCO\n+o4D1O1plyspN4FUPLsaaPnOLf8xncitPsbqNFHmsVHbWkqZx4arxkpptXVDj+DarmTgbXP5+bp/\nYmDwj7BYvOzr+AzFxf41+dqZTIYXX3yRl156Cb1ez7lz5zh69CgGQ2G+rTJTMWIXpkl0ByGnYfI7\ncJ6rx9LmRjFun25OU1VmBvsZvdbF6LXLzI8OA1BkteFt3099xwF8HQdwlMt5uELJZVQiswkWpmL5\n4ciVj/HFV/exmSwGyjxWmjorKauxUuax4qqxbdo9bNuBDLxtTNMy9A/8V6anv0xZ2Rn2tv3lmmw5\nEEJw69YtnnrqKaLRKO3t7Zw7d64gJ6QIVSN5Y4HYhWky40soRh3WQxXYTtRsq83h8cUIY91XGL3a\nxdj1K6TjcRSdjpqW3Zz82Q9R37GfygY5TLnRhBAshVIsTMYITcVWO7fofILbh1jpDTpKq4upbS3F\nVWOlzGOjzGPFWmKW86YbTAbeNpXJhLne8wmi0S589Z+goeG3UZR7fzFcWFjgySefZHh4mIqKCj7y\nkY/g8/nuveC7pC5niF+aJXZxBm0pg95VhPOJBqydlei2wSo0TVOZHRrMd3FXLzM3MgiAtaSUpsPH\n8e/vpL5jP0VWuR9uo2Qzar5jm4wRmoythlw2tbJJWwFnuYUyj42mzgrKavLB5iy3oNNvnxGGrWzr\nvzJIPyKRGOVa98dIp2doa/u/qKr8qXv+mrlcjgsXLnD+/HkMBgOPPvoohw8fRr/BV7hkJpbzw5bX\ng6AKzM0l2N7ZRFGrC2WNz+jbaIml6Gu6uKuklpdQFB3Vza2c/NkP4T/QSUW9X97qvc6EEMQi6dVg\nC03mO7fF+UT+mmryp4+4PTZaj1bhrrVRVmujrMYmV0ZucjLwtpnFxS6u9/wKAAcO/OOa3F0XCAT4\n5je/STAYZM+ePTz22GPY7Rt3VqLQBKm+BZZfmCIztoRi0mM9UoXteM2WvlBVCEFwfJSRy5cYufIK\nM8MDIAQWh5OGA5349x+ift9BeS7lOlKzGuGZOKHJ5dd1ben4q4tIHO4iyjw2mjsrcNfaKau14Sgr\n2vJvsHYiGXjbyOzcN+nt/T0sFg/7Oj5LcbHvnr5eKpXimWeeoaurC4fDwfve9z5aW1vXptg7oGVU\nElfmiL04TS6URF9ixvm2lWHLoq35ravmskzc7GH48iWGL19kORQERaGqsZkT73k//gOdVPobZRe3\nDlLxLMGJZYKBZUITsdWl/2JlT5vBpKPMY6PxYAVujy3fuXlsO3qj9nYj/09uA0IIxsf/J8Mjf0aJ\n8zAdHX9zz5vJe3t7efLJJ4nH4xw7dowzZ85gNpvXqOKfTF3OEHt5mvgPZtASOYy1Nlzv35VfbbkF\nD29OLi8xeu0yw10XGeu+TCaZxGAyU99xgOPvfh8NBw/Lo7vWWDyaXgm2ZYKBGMHAMsvh1Oqf20rN\nuGtt+Pe5cdfacdfacJRb1vzqGmlzkYG3xQmhMTD4h0xOfpHKyp9mz+4/Rad768EUj8f59re/TW9v\nL5WVlbzvfe/D49mYG4mz8wliL0wRvzoHqqBolwv76VpMfseWW80WmZliuOsiw5cvMXWrFyE0rKUu\ndp24n4ZDR/C278No2pg3ENuZEILlhRTBwPJK9xYjNLFMYimz+piSymKqGhzsvd9DuddOeZ2dIptc\n+r8TycDbwjQtQ2/f7zE390283l+kqfFT9xQMvb29fOtb3yKVSvHggw9y8uTJDVmUkplYZunZCVK9\nC2DQYT1Uie2UZ0ttEheaxuzwIIOXLjDcdZHw9CQA5fV+jr7zvTQeOkplQ5McqrwHmiZYnEusdG35\ngAtNxFY3bSs6BVe1Fe8eF26vnXJvvnMzbdHhb2ntye+ELUpVE1zv+VXC4Rdoavw96ut/+S1/rXg8\nzpNPPsnNmzeprq7mwx/+MJWVlWtY7Y8SQpAejrL83ATpoUUUiwH7WS+249Vb5joeTVWZ7LvJ4KUL\nDL3yMrHwAjq9gbq2dvY/8gSNh47Kzd9vkaYJIjNx5sdXwi2wTGhymVwmf9SW3qCjrNZG06GKfLDV\n2SnzWDFs88O/pXsjA28LymYXudb9iywtdbN7159QU/Mzb/lr9fX18a1vfYtkMsmZM2c4derUunZ1\nQghSfWGWn5sgE1hGZzfifNyP9WgVOvPm/3bMZbMEblxj8OLLDHf9gOTyEgaTGd++g7S8/wT+g4fl\n3ri7JDRBNJhkfnyJ+fFl5seXCAZeDTdjkR53rY09p2pWhyRLqorRy71t0l3a/K8w0uuk0/Ncvfbz\nJBLjtLf/31SUP/KWvk4qleLJJ5/k+vXrVFVV8aEPfYiqqqo1rvZVQhMkrwdZfm6C7GwCfamZknc0\nYj1UtemP/cqmUox2X2bw4gVGrlwik0xishTTcPAwLUdP4tt3EGNRUaHL3BKEECyHU8yPLRMMLDE3\nlu/eMsn8sKTBqMNdZ2fPyRoqfA4q6u2UVBTLLQDSmpCBt4Wk0rNcvfpB0uk59u/7HC7Xibf0dQKB\nAF/72teIRqPcf//93HfffevW1QlNkLg2z/L3J8iFkhjKLZS+t4Xi/eUom/gdeiaZYPjyJQZ+8CJj\n166Qy2aw2B20HDtNy9ET1O3dh6HAd/ttBfFomvmxVzu3+fFlUrEsADq9grvWRvPhSirq7VTUO3BV\nF8tTSaR1IwNvi0ilprly9QNkMmH27/s7Sko67/prqKrK888/z/PPP4/T6eSjH/0oXq93Hap9taNb\n+l6AXDCJsdqK6wO789fybNJ369l0ipErXfS//DyjV7rIZTPYSl3sffBhWo6ewLOrDd0GnyyzlWSS\nOebGl5gbXcqH3NgS8Wh+taSigKvGiq/DTWW9nQqfg7IaG/pN3t1L28ubBp6iKH8mhPgPG1GM9OMl\nk5NcufpBstkIB/Z/AafzwF1/jXA4zNe+9jUmJyfZt28fjz32GEXrMAwnNEHyRoilZ8bJzScxVBZT\n9sHdFO3ZnEGXy2YZu3aZ/pdfYLjrItl0imJnCXsffJjWE6fxtOyWKyt/jNuLSmZHosyN5UMuPBNf\nPXqrpLIYT2spFfX5YUl3nV0euyUV3J10eA+uexXSG0omA1y58gFyaoyDB76Iw9Fx11/j+vXrfOtb\n30JRFN797nfT3t6+5nUKTZC8uZAPurkEhori/Gbxve5NF3RqLkeg5xq3LjzP0Cs/IJNMUGR3sPvU\nA7SeOE3tnr3y1oEfEo+mmRtdWgm3KPNjy2TT+UOTzVYDlT4nTYcqqPQ5qPA55BU30qYkhzQ3sVRq\neiXsEhw48EUc9r139fnZbJannnqKy5cv4/V6ede73kVJScma1iiEINUfYem7Y2Rn4hjKLbh+rhVL\nR/mmCjqhaUzeuknfi88xePECqdgyZquV5qMn2HX8NHV796Ev0D1+m00uqxKaiDE3usTsaJS50SWW\nF/KnlOh0Cu46G7uOVVHZ4KTS58BZYdlyBwNIG0eNRslMTpKdnCI7OUl2anL19+5PfALnT71tw2q5\nk5/wfYqijAI9wI3XfLwlhMiuZ3E7WTo9z5WrHySnLnPgwD/eddgtLCzwL//yL8zOznLq1CnOnDmz\n5gtT0oElot8ZIzMaRe8qovRnWynet7mCbmEyQO/z36fvxfMsLwQxFllo6jxK64nT1HcclAtPgFgk\nxcxwlNnhKLOjS4QmltHU/NikzWWmyu+k40wtlX4n5XU2DPLmbek1tGSS7NTUjwm1/D9ry8uve7zO\n4cBY68Hc2Ih+jd+Av5k7CbzrwBPAXqAdeBj434BmRVEmhBB390osvalMJszVaz9PJhPkwP6/v+uw\n6+3t5etf/zqKovD+97+flpaWNa0vG0yw9NQYyZsL6GxGSt7eiPVwFYphc8x1xSJhbr10nr4XnmN+\nbBhFp8O37yCnP/ARmjqPYjTv3C0EmiZYmIoxOxxlZjjKzPAisXD+Fm6DSUdFvYP9D9VR6XdS6Xdg\ndcrjz3Y6oWnk5ubITEyQnZhY+ThJdnKSzNQUaij0uscrRUUYPR6MtR6KDxzAWFuLsdaDqbYWY20t\n+gJcFH3bHY3hCCGmgWng32//OyU/htG0TnXtWLncMte6P0IyGWDfvs/hdB68i8/N8fTTT3Px4kU8\nHg/vfe9713QIU11Ks/RMgHjXLIpBj+MhL7bTteg2wWKETCrJ0KWX6X3hWQI93QihUdXYzJmPfJxd\nJ+6j2Lmx7yQ3i0wqx9zIEjPDi8wM54cnb8+9WZ0mqhpL2H/WSXWTk7Jam9zMvUNpyWQ+yCYnyQQC\nZCcmyUwEVoNNZF8zmKfXY6ypwVjrwX7mAYyefJCZaj35QCsr27RD3HcSeH/94/6lEEIAg2tbzs6m\nqkmudX+MWGyAjva/wVV6/I4/NxaL8ZWvfIVAIMCxY8d46KGHMKzRnJSWVlk+P0HshSmEJrAdq8H+\nYF3BjwATmkbg5nVuPvcMg6+8TC6dxlFeydF3vpddpx6gzFNX0Po22u1N3a92b1HCUzGEABQo89ho\nPVZFdaOTqkYndlfRpn1hktaWEAI1FHq1SwtMkJ3Mf8xMTqAGX9+l6Ww2jN46zC0t2M8+iLHOi6mu\nFmNdHcbqapQtOt/9plULIT67EYXsdJqW48bN3yIavcLevf8Dt/vMHX/u1NQUX/7yl0kkEmu6ClNo\ngsTVeaJPjaEtZ7DsK8f5cD2GMsuafP23aik4z43nnuHm+WdYCs5jtlrZc/oMu0+fwdO6Z8e8iAsh\niMwkmB6MMD24yPRQlPjiyvCkWU+V38Ghx335gPM75b1u25yWyeTn0CYCZCYmf+SjSL16PRKKgqGq\nClNdHbb77sNU58VYV4vJ6813aSUl2/LnSP4EbAJCCAYG/4BQ6BlaWj5NZcVjd/y5169f5xvf+AZW\nq5WPfexjVFdXr0lN6bEoi98aITsZw1hnp+yDuzHXF27sPZfJMPTKy9x47hnGe64BUN++n9Pv+zBN\nh49jMG2NA6fvhaYJFiZj+XAbXGR6aHH11JJip4ma5hKqG0uobnRS5rHKE0u2IS2ZzHdlgXGy4+Nk\nxgNkxsfJBALk5ubIt/N5isWSnzer82I9cRJjXR0mbx3G2jqMtR50O+Bn5ofJwNsExsf/hqmpL1Hv\n/WXqan/+jj5HVVWeeeYZXn75Zerr6/mZn/kZrFbrPdeSi6SIPjVGsjuI3mEq6MpLIQTzo8PceO5p\n+l58jnQ8jqO8khPveT9t95/d9jcRqKpGMLDM9EA+3GaGoqtnTtrLivDtLaO6uQRPSwkOt9wasF2s\nhtr4GNnASqCtBFtubu51j9WXlmKqr6f4yGFMdd58oNXVYaqrQ+92y++JHyIDr8BmZv6V4ZE/p6ry\n7TQ23tmBNqlUiq9+9asMDQ1x5MgRHnnkkXveciCyKkvPTbJ8Pn+Pm/2sF/v9tegKsAQ9nYjT+8Kz\n9HzvuwTHR9EbjTQfOcHeM+fwtnVs25NPclmV+bElpgcXmRpYZHYkunpjQEllMU2HKqhpLqGmuQS7\na+euNN0O8qF2O8zG88E29ppO7TX0LhcmrxfrsWMY672Y6usxeesx1XsLuuJxK5KBV0Dh8Ev03fpP\nlJaeYPfuP0VR3vyFPBqN8qUvfYn5+Xne9ra30dl592dq/rBkf5jFrw+jhlP5ebrHfBhKNv4FdW5k\niO5nvsOtF8+TTaeo8Ddy9mO/yq4T91Fk235X7qiqRnB8mclbESb7I8wOR1Fz+YAr89jYfaJmNeCK\nHTtv+Gmre8NQGx8nNz//usfqXS5M9fVYjx3D5KvH6PViqvflQ81uL9DfYPuRgVcgicQoPTd+neJi\nPx3t/w863Zu/oM3MzPClL32JdDrNBz7wAZqa7m1XSG4xTfSbwyRvLmAot+D+xXaKmjZ2+X42neLW\nhee5/vR3mB0exGAys+vk/ew79xhVjc0bWst6E5ogNBVjqj/C5K38QpPbWwTKam3svc+Dp7WE/5+9\n9w6Po7oX99/tfSWtepesZlmyLDfcCwZsMGBMDS0EAiEhIYEn+abdm/vLDWn33iQ3JLlpBAIJJSFg\nmk3HYGPjIluusopl9V5WK2l7nd8fs1pJWDauKta8z7PPzM6enZk1q335nPM5n5OcGy2V5poiCKEQ\n/vYOfA0N4qOxEV9jA96GRgIdHaPaKmJjxUht6VLU4UhNNRSpSVIbFyThTQCBgJ3DR76MTCZnTskT\nKJWf/WWvra3lpZdeQqvV8sUvfvG81q4TgiEcO9sZ3NoEApjXZWFakTquE8etrc0c/uBtKrd/iNfl\nJDYtg8vv/TKzVl5+ySygKggC/V2uiODajvfjcYpJJtGJegoWJZFaEENqQTS6KbLK+3QlODCAr6EB\nb2MjvobGYcE1NSH4fJF2cqMRdXY2+oULUGdlocnKQpWZiTpDktpkQBLeOCMIQSqOPYLb3cTc0r+h\n03328jzl5eVs2bKFxMRE7rzzTszn0W/vbRzA9uoJAl0utIUWoq/PQTlO40GhYJAT+/dw8J3NtFZW\nIFcoyV+8jDlXXUPqzKJLYoDd3ucJd1H20VZtiyyPY4zRkFUSS1pBDKkFMRhjpDG4yYbg8+FrbY3I\nzNvQIHZBNjQQ7OsbbqhUok5LQ52djWH5ctTZWWiys1FnZ0/qSdcSkvDGnRMn/hurdTsFBT8mJmbx\nZwEmzrUAACAASURBVLbfsWMHW7duJTc3l1tvvRWN5txKPYW8QQbeacC5pwNFlIbYe2ahmxV7Tuc6\nWzwOB0c/eo9D725hsKcbc3wCK+68l+LVV075Cih+b5C24zZaqvpoqezD1ukCQGdSkVoQExFcVLyU\nRTlZCPT14aurw1vfMCy3xgb8rW0QDEbaKeLiUGdlYrpiDeosUWjqrCzU6WnIpBqsUxJJeONIZ+cb\nNLc8RVrq50lLvfO0bQVB4P3332fXrl3Mnj2bjRs3nnMmpue4DdsrtQQHvBiXpGBelzUu5cCsbS0c\nfHszxz7eSsDrJW1WMavveYCcBYum7PI7Q+NwLZV9NFf20VHXTyggoFDJSc2LZtbyFNILLVhSDJLg\nJhBBEAh0dOCtq8dXX4e3rh5vXR2+ujqC/f2RdjKtFnVmJtrCWZjXr0eTlRURm5QBeekhCW+ccDpP\nUF3z70RFzScv799P2zYUCrF582YOHjzIwoULueaaa5CfQyp+yOWn/80GXOVdKON1xH+5BE1W1Ll+\nhDNCCIVoPHyAA2+/QePhAyhUKmYuW8W8azaQkDXjol77YuEc8NJa1UdzOIpz28VxuNhUAyWXp5NR\naCE5N0paRWACEAIBcc5afR3eE3V46+vw1dXjbWhAcLki7RTR0ahzcjBddRXqnBlocnLRzMhGmZx8\nyU5zkTgZSXjjQDDo4mjFw8jlWoqLf4tcfurukEAgwKZNm6iqqmLlypVcfvnl5xQpuCut2F6tJeT0\nY1qdjvmKDGSqi/eHHfD7qfz4Q/ZveRVbeyuG6BiW3nYXc668Zsp1WwaDITpPDNBUYaW5qg9rqwMQ\nuynTCy2kz7KQXmiRVhIYR0Iej9j1OFJq9XX4mpphRGFjZVISmhkziL75ZjQ5OWhyZqDOyUFpsUzg\n3UtMFiThXWQEQaC6+j9wOk8wt/RvaDWnzq4MBAL861//4vjx46xbt44lS868ePQQIW+A/s31uPZ3\noUo2EHdvMerUi5f16HU5Ofz+2xx463Wc/TYSsnJY//C3yF+yHIVy6oxzuAZ9NFVYaaqw0lLVh88d\nQK6QkZwbxZIbc0gvtBCXZpxUa/1dioQ8HrHr8cQJvLW1eGpr8dXV429rGy6bJZejTk8XI7bLL0ed\nk4MmJwd19gwUxvOvNiRx6SIJ7yLT3v5POrteIzv7USyWZadsFwgEePHFF6mtreXaa69l4cKFZ30t\nb+MAff86TtDmwbQ6DfOVmRdtqoGjz0r5W69z5IO38bndZMwu5ZqvfYuM2XOmxNiVEBLobrbTdLSX\npgor3U3iIpWGKDW58+LJnB1H2swY1FrpT+RiIPh8Yor/iRN4amvx1tbiqz2Br6UFQuLke5lKhTo7\nG13JbKI2bkSTm4N6xgzUWVnTsg6kxPkj/TVfRJzOExyv/TGWmOVkZ33tlO1Gyu5cqqcIgRCDHzRh\n396KIkZ7UcfqrG0t7N/8KlU7PiQUDJG/eBkLN9xM4ozJvzSi1x2gpbKPpopemo714R70gQySss0s\n2pBNZnEccenGKSHsqYIQDOJvaRmWWjhy8zY0QkCsC4pCgTozE01BAebrrkOTl4cmPw91RsaUXYZG\nYnIifZsuEqGQj2PHvolCYWDWrF+csmzY+crO3+2i7x/V+Duc6BckEn39DOSaC/+ftae5kT2b/snx\nvZ+gVKqYfcU65l97I9GJ5z4Bfjyw93loONxD/aFeOmr7CYUENHolGbMsZM6OI6PIIk36vgAIgkCg\nvV3sgvxUd6Tg9UbaqdLT0eTmYrx8jSi2vFzU2dnIz3G6jYTE2SAJ7yJRX/9r7I5jlMz+ExrN2FX9\ng8EgL7300jnLznmgi/7XTiBTyon9/Cx0RRd+Xl13Yz17Nv2T2rJdqHU6Fm28lXnrb0BvvrjZnueK\nIAhY25xhyfXQ2yImnMQk6Sm9Kp3M4jiSZpilpXPOg5DLhffECTzV1XhrjuOpEbchuz3SRpmYiCYv\nD8OixRGxaWbMQH4BVvSQkDhXJOFdBGy2PTQ1/4WUlM8RH3/VmG1CoRCvv/46NTU1rF+//qxkF/IF\n6X+9Dld5F+osM7F3zERxgTMGu+pPsHvTP6nbvweN3sDim+9g3voN6IyTrzxSKBiio26AhkO9NBzp\nYbDXE+6qjGLJTTnMmBNPdKJ+om9zyhGJ2mpq8NbU4KkWt76mpkgCiVyvD3dFXou2oABNfj6a3Fxp\nDpvEpEQS3gUmELBTWfltdLpM8vN+MGYbQRB45513OHLkCGvWrOGyyy474/P7u5xYn68m0OPCtCYd\n8xWZyBQXbsyps66W3S+/QP2BfWgNRpbedhdzr75+0tW3DPiCNFf20XC4h8YjVjxOPwqlnLTCGOZf\nnUVWSZy0wsBZEHK5xG7Imhq81TV4jtecFLWpMjLQFuRjvu46tDML0BQUoEpNleaxSUwZJlx4MplM\nAewH2gRBuG6i7+d8OVH3SzzeDhbMfwmFYuyoYtu2bZSVlbFkyRJWrFhxxud2loe7MDUK4r5YjDYv\n5kLdNtbWFj558Vlqy3ahNZpYfvs9lK67Do1+8kRGfl+Q5gorJw5003jUSsAbRKNXkjk7lhlz4kmf\nZZGyKs+AQE8PnqoqPJWVeKqqPztqKyhAk5cvpfxLTHkmw6/DI0AVMOX7QGy2MtraniM9/T6iouaO\n2aasrIzt27dTWlrK2rVrzygjUAiGGHizAceudjQzorDcMROF6cJEL4M93ex66QUqP/4QlVbDklvu\nZP61GyeN6PzeIE0VVk6Ud9NU0UvAF0JnUlFwWSI58xJIyY9GIY3HjYkgCPjb2vFUHosIzltZRaCn\nJ9JGlZ6OdmY4O7IgH+3MmVLUJnHJMqHCk8lkacC1wE+Bb07kvZwvwaCHqurvo9WmkzNj7I9SU1PD\n22+/TUFBAddff/0ZyS5o92F9oQpfwyDG5alEXZN9Qbownf029r72L468/zbIZMxbv4HLNt46KZJR\nfJ4ATRVW6g5001RhHZbc4mRy58WTkhctJZ18CiEYxNfUhOdY5YjorYrQwIDYQC5Hk5ODYekStLNm\noSksRFtYKC1ZIzGtmOgI73HgO8CU/6traPgtbncjc0v/PmZXZnt7Oy+//DLJycncfPPNZ1QI2tdi\nx/pcJSFXAMvnCtDPHTvb82zwedzse2MT5VteI+D3UXz5VSy5+Q5MsXHnfe7zIeAXI7nasi5Rcv4Q\nOrOamYuTyZmfIEpOqnIChCdt19WJUhsSXHU1gtsNiBO2NQUFmNetQzurUBRcfj5yrbQkkcT44w64\n6XZ10+3qpsfVI+67xec35t7IstRTF+S40EyY8GQy2XVAtyAI5TKZbPVp2j0IPAiQkfHZa8dNBHZ7\nFc0tT5KSfNuY1VT6+/t54YUX0Ov13HHHHajPoEqE82A3tk3HURjVxH9lznmXBwuFglRu/5CdLz6L\n09ZH/pIVLLvtbiwpqed13vO7J4H24zaOl3VRd7AHnzuAzqRi5tJkcuclkCxJTozc6utxVxzDc/Qo\n7ooKvNXVkUVH5Xo9msJCom+5BW1hIdqiWWhmzJCWr5G46PhDfqxua0RmEam5e0YJzu63n/RenVJH\ngj6Bfm//GGe+eMiEofp044xMJvs58HkgAGgRx/BeEQTh7lO9Z8GCBcL+/fvH6Q7PDEEQKD9wOy5X\nPUsWf4BKNbpL0Ov18tRTTzEwMMD9999PQsLpozRBEBj8oBn71mZxvO6uQhSG8/vxaq44wrZnn6Sn\nsZ7kvAJW3/MlUvJnntc5zxVBEOhtcVBT1smJfV04B3yoNApmzI0nf2EiaTNjpm13pSAI+FtbRbEd\nrcBz9CieykpC4ar/cr0ebVER2tmz0RbNQjtrFurMTGm8TeKCEhJC9Hv7TxmVDT3v8/QhMNofSpmS\neH088fp4EvWJxOtG7OvjSdAlkKBPwKC6cMtnyWSyckEQzmhe14RFeIIgfB/4PkA4wvt/p5PdZKWr\n6w0GBvZTOPPnJ8lOEARee+01enp6uPvuuz9bdoEQtk21uA52o5+fSMyNuedVC9PW0cb2556mbv8e\nTHHxXPuNb1OwdOWElM4a6HFTu6+T42Vd2DpdyBUyMopiWXZZIlklcaim4dI6/q5uPBVHcR89iqfi\nGJ6KishabTK1Gk3hTKJuvBHt7GJ0xcWos7ORneOaiBISAMFQkF53L12uLvHh/NQ2/AiEAie916K1\nRMRVFFdEgi5BlJhelFi8Lp4YbQzyU1SVmgxM9BjelCYQsFN74ueYTSUkJ99y0us7d+6kqqqKtWvX\nkpOTc9pzhVx+ep+txNcwiHltJqbL089ZTH6Ph92v/JPyLa+hVKtYfscXmLd+Ayr1+JZv8nkC1B3o\noXp3B+214g95Sl40c65IJ2deAtrzjFynEiGnE3fFMdyHDuE+cgTP0aMEurvFFxUKNHl5mK66Em1R\nMdrZxWjz8pBJBZIlzgJ/0E+3u3tsiTm76HR1YnVbCQrBUe/TKDQk6hNJNCQyN2EuCfoEEvWJosTC\n0VmsLhbVaZY1mypMCuEJgrAN2DbBt3HWNDT+Hz5fL3NKnjipVmZtbS1bt26luLj4M5f5CfR56H26\ngkCfB8vtBehLzy05RRAEast2se1vT2K39lC0+kpW3PEFDNEXbr7emdxDx4l+qnZ1cOJADwFvkKgE\nHYs3ziD/siRMlks/cUIQBHwNjbgPH8Z9+BDuQ4fxHj8eWQVAnZWFftEidLOL0RbPRls4E7lON8F3\nLTGZGUr8GCmxTmdnRGbdrm6sHutJ79Mr9SQZkkjUJ7I0ZWlEbIl68ZFkSMKsNk+bgumTQnhTEZer\nkZaWZ0hJvhWzuWTUa1arlU2bNpGYmMiGDRtO+2XydznpfaqCkC9E/AOz0WSf27QAW0cbHz79ZxoP\nHyA+M5trv/FtUmfOOqdznQv2Pg81ezqo2t3JYI8blUZB3oIECpckk5QTdUn/QQXtdtxHjojR2+HD\nuA8fiUwHkBuN6EpKMH3lK+hK56ArKUERPbUWxJW4uPiDfjpdnXQ6O+lwdkS2I+U24B046X1RmqiI\nuIriiiL7I4VmVE+uCkkTjSS8c6Su/lfIZCpmfGrOXSAQ4KWXXgLg9ttvP21Gpq/FTu/TFaCQkfCV\nElRJZ1/Jwu/zUvbqv9j3xiYUKjWX3/sgpWuvRT4OYz3BYIiGQ71U7myjpdoGAqQWRLPw2ixy5iag\n0lx6401CKISvvh7XwYOi3A4dwldXL1YpkcnQ5OZgXnsVujlz0M2ZgzonR0oqmcYIgkCfp+8kmY3c\nt7qtJyV/WLQWkgxJpBpTmZcwb5TEEg1id6NOKfUKnC2S8M6BwcEjdHe/RVbWw2g08aNee//99+ns\n7OT2228nJubUXYmeEzasf69EblQTf38xytiz//K2Vh/jvT/9FltHG4XLV7Py7i9ijLGc9XnOlsFe\nN5U726na1YFr0IfRomHh+ixmLknGHHdp/REKPh/uY8dwHziAa3857gMHCA5Fb1FR6OaUYF6/XhRc\nSYk0kXua4fK7xOjM0UmnKywzR8cowflCvlHv0Sl1JBmSSNInsTJtJUmGJJINySQbkiPdj1rlpd/1\nPxFIwjtLBEGg9sR/oVJZyMz40qjXampq2Lt3L4sWLWLmzFOn/bsrrVifr0IZpyP+/mIU5rNLJvG5\nXez4x9849O6bmOMTueXff0JmSek5fZ4zJRQSaKqwcuzjNpqOWZEBmbPjKFqRQkZR7CUzXy5ot+M+\ndAhXeTnu/eW4jx6NrOemzszEeMUV6OfPRzd3LursrEu6q3a6IwgCVo+VNkcbHY6OSGQ2JLJOZ+dJ\n88jkMjnxuniSDEnMip3FFRlXiHIbIbUozaXdxT+ZkYR3llj7ttPfv5f8/B+iVA73jw8ODvLaa6+R\nlJTEVVeNvSQQgLtKlJ0q2UD8F4uR688u86nhUDnv/+X/sFt7mXfNBpbd/nnU2osXVTn7vVR+0k7l\nznYcNi/6KDULrsli1vKUSyIBxd/Vjbt8P67yA7gOHMBbUyMmlygUaAsLibn9dnTz56GfNw9l3MRW\no5G4sISEEL3uXtod7eLD2U6boy3yvMPZgTfoHfUek8pEklGU15z4OSfJLF4ff0lkM16qSMI7CwRB\noKH+N2i1aaSm3B45HgqFePXVVwkEAtxyyy0olWP/s7qr+7A+F5bd/bOR6878n9/rcvLRM09wbPtW\nLKnp3PHY/5CSX3jen+lUdDUMcvjDFurKuwmFBNJnWVhxWz6ZJbFTulizv7MT1969OMvKcJXtw9/S\nAoBMp0NXOoe4r34V/fx56EpKpMVKpzghIUSPq+ckkQ3Jrd3Rjj/kH/Uei9ZCiiGF/Jh8VqevJsWY\nQqoxNSI0KQlkaiMJ7yzo6/uYQfsRZhb8FLl8OBll//79NDQ0cP311xN3iijAXd2H9dnKc5Jda2UF\nb/3+Vzj6rCy68XMsvvl2lBehdFQwGKL+QA+HP2yhq2EQtVbB7NVpFK9OJTphcqyecLb4u7pwlZXh\nKivDubcMf3MzII6/6RcuIOauO9HPX4B2ZoFUjmuKEQwF6XH30O4YIbQRcutwdpw0gdqitZBqTGWm\nZSZrMtaQakglxZhCijGFZEMyetXU/J5LnBmS8M4QQRCob/gdWk0Kyck3RY739fXx/vvvk5uby7x5\n88Z8r6fWJsouKdyNeYayC/j97PrXc+zb/ArRiUnc8dgvSM4ruCCfZ9T9Ofwc29nG0W1tOPu9RMXr\nWPG5PGYuSZ5y68v5u7txle0TJbd3r7jOGyA3mdAvXIjlrjvRX3YZmoICKXtyCuDwOWhztNFib6HV\n3kqrozWybXO0nSS0OF0cKcYUimKLuCrzKlKNo4UmZTZOb6bWr9kEYrPtYnDwIAUFP45Ed6FQiNdf\nfx25XH7K5X58LXZRdvF64u8/8zG73pYm3vrdL+lpaqDkiqtZdc/9F3ysbqDHxcH3W6je3UHQHyJt\nZgyr7yogsygW2RRJQgkODODcsxfn7l249pbha2gAxPlv+gULiL79dvSXLUQ7c6ZUlmsSEgwF6XJ1\njZbZiH2b1zaqvVltJs2URkFMAVdkXEGqMTUitWRDspTdKHFaJOGdIQ0Nv0OjSSIl+ebIsX379tHU\n1MQNN9xAVNTJE8b93S56n65AblQTd4ayEwSBox++x0dP/xm1Xs/G7/wHOfMXXdDP0tNs58C7TdQd\n6EamkDFzURIlV6QTmzL5xycEnw/34cM4du3CuWsXnqMVEAoh1+vRLVxA9C03o79sEdpZhZLgJgl2\nn/2UQmt3to+K0pQyJcnGZNKMaVyZeSVppjTSjGmkmdJINaYSpZn49Rolpi6S8M6AgYFD9A/sIy/v\nB8jl4hSCwcFBtm7dSk5ODqWlJ08JCPR76X3qKMhl4tSDM1ih3Odx88GTf6Bqx0dkzC5l/cPfumBl\nwQRBoLXGxsF3m2ipsqHWKpi7NoOSNekYosa3xubZIAgCvvp6nJ98gvOTXTj37UNwuUAuR1dSQtxX\nvoJh2VJ0JSXSGNwEIQgC/d5+mu3NNA8202xvpmmwiZbBFlocLSdVCYnSRJFmTKMwtpCrMq8SpRYW\nW5IhCaVc+lmSuDhI36wzoLnlryiVJlKSb40ce/fddwkGg1x77bUndWWG3AF6/1pByBsk/sGSM5pU\n3tvcyOZf/xd9HW0svfUuFt10G3L5+UcoQkig4XAv5e800t1kR29Ws+TGHIpWpqI5i8SZ8SRgs4ly\n++QTnLt2EejqAkCVmUHUDRswLF2KYdEiFGbzBN/p9EEQBGxeW0RozYPDcmu2N2P3Da95JkNGijGF\ndFM66zLXjRJaqikVs1r67yYxMUzOX7xJhMfTTk/PO6Sn3xeZd1dXV8exY8dYvXo1FsvoyiZCUMD6\nQhUBq5u4LxajPoNuwsodH/H+E/+HWqfj1h/8hIziOed930OiK3uzAWurg6h4HavvKqBgcRJK1eTq\n6hMEAU9lJc6PP8ax/WPcR46I3ZRRURgWL8awbCmGpctQp03cYrXTgU9LbShKa7KL25ELecplcpIN\nyWSYMlifvZ4MUwaZ5kzSzemkGdNQK6SVHiQmH5LwPoOW1r8BkJ72BUCslfnmm29isVhYtmz06uaC\nINC/uQ5vbT8xt+ShzTl9keBQMMjHLzxD+ZZXSZtVzHWPfPe8uzAFQRTdvjcb6G1xEJWg48r7ZpG3\nIGFSLawatNtx7tqN4+PtOD/eQaCnBwDt7NnEPfQQxpUr0BYXS+NwFwGX30XjYCONA42RbZO9iebB\nZhx+R6SdXCYnxZBChjmDkhklZJjDUjOlk2pMlaQmMeWQhHcaAgEn7e0vEh9/NVptCgC7d++mr6+P\nu+++G9Wnxoycu9px7unAuCoNw4Kk057b43Cw5Tf/TdORg5Suu47V9zyA4hQT1s+EMUV3byF5CxMn\nhegEQcBXV4dj+8c4Pv4YV3k5BALITSYMy5dhXLkK44rlUjWTC0QwFKTD2XGS2BoGG+h2dUfaDXU/\nZpmzKJlRQqY5kwxzBhmmDFKNqagU0rioxKWDJLzT0N39JoGAnfS0ewBwOp3s2LGDgoICcnNzR7X1\n1PXTv6Ue7axYotZlnfa81tZmXvvFjxns6eGqB79OyRXrzus+22ps7HrlBN1NdqLidVxxbyH5k0B0\nQjCI+9Ah7B9sxf7hVvxN4qRvTX4+sffdi3HVKnSlpcjOQ/TTHbvPHhFaw0BDZNs82DyqaLFJbSLb\nnM3i5MVkmbPIisoiy5xFhjkDjWLyJi1JSFxIpF+a09DW/iJ6fS5RUfMB2L59O36/nyuvvHJUu+CA\nl75/VKOM12H5XMFp57A1VxzmjV/9DIVKxW0//DmpBedeHqy31cHuV+toPmbFGKNhzT0zKViUNKGi\nC3k8OHftxv7hVhwffkSwrw9UKgyLFxN7330YV61ClZw8Yfc3FREEgS5XF/X99dQN1FE/UC/KbaBx\n1KKfCpmCNFMaWeYslqUsIysqi+yobLLMWVi0FqlgscS0RxLeKXA4ahgcPERe7r8jk8no7e1l//79\nzJ8/n/j44SWBhGAI6wvVCL4gsQ+WID/NGnCVOz7i3T/+hpjkFG76/n9ijju3lc3tfR7K3qinem8n\nGp2SJTflULI6DaV6Ysa7gv39OLZvx/7BVhw7dyK43ciNRoyrVmG6Yg2GlStRGCf/HL+JJiSE6HB2\nUNdfNyy38Nbpd0baRWuiyY7KZmXaykiklhWVRboxXeqClJA4DZLwTkFb+z+RydQkJW0EYOvWrSgU\nClavXj2q3cDbjfiaBrHcUYDqFPUmBUGg7LWX2PnPv5M2q5gbvvUDtOcgAJ8nQPnbTRzeKhY8nntl\nBvOuzkRrGP8fuYDNhmPrVgbffgfnnj0QDKJMTCT6xo0Y11yB4bKFyE6z+O10JhgK0upoFcU2UE9d\nfx11/XU0DjbiDrgj7eJ0ceRE5bAhZwM5UTnMiJ5BTnQOFu3FX/NQQuJSRBLeGASDHjo7XyMhfi1q\ntYWOjg6qqqpYvXo1xhGich/rxbGzDePSFPRzxo7WhFCIrU//mcPvvcnMZatY99CjZ134WRAEjpd1\nsfuVEzgHfOQvSmTxDTnjvjxPsL8f+0jJBQKoMjKI/eIXMa29Cm1RkVSfcgQhIUSbo41aWy21tlpR\nbAN1NA40jhpfS9QnkhOdw/zE+eRE55ATncOMqBlSVREJiQuMJLwx6O39gEBgkJSU2wBx7E6r1bJ4\n8eJIm6Ddh+2VWlQpBqLWZ495nlAwyLt/fJzKHR8x/7obWXXXfWcthO6mQXa8WEtn/QAJmSau/vJs\nkmaM3w9hcGAA+wdbGXznHZy7d4uSS08n9r77MF29Du2sWdLYEGDz2ESx9ddGBHei/wSugCvSJtWY\nSk50DstSlpEdlR0Rm7TkjITE+CAJbww6uzajVicQE7OYzs5OqqurWb16NVqtGFEJgoDt5eOEvCHi\nP1eATHmyxIIBP2/99pcc3/sJS2+7i8U33X5WYvA4/ex+tY7KT9rRGVVc/vmZFC5JHpeiziGPB8e2\nbQy8/gaOnTvB70eVmkrsvV/AdPU1aIumr+Q8AQ91A3URqQ1JrtfdG2kTrYkmLyaPjbkbyYvJIy8m\nj9zoXAwqaX09iWlIwAuuPnD3gcsq7rus4vOcKyB17FVmLgaS8D6F3z+I1foxaWl3IZMp2L59OxqN\nhkWLhgs4O/d24KmxEX39DFSJJ/+IBXw+Nv/659Qf2Meqz9/PgutuPOPrC4JA7b4udr5Ui8cZYM6a\ndBZem4XmLFdGP1uEUAjXvv0MbH4D+zvvEnI4UCYkYLn7bszrrxEngU8jyQ1lRtb01VDVV8Vx23Fq\nbbU025sJCSEANAoNM6JmsDRlKfkx+eRFi3KL08VNq38riWmEzzWGuGzDz0e9Fn7uc5z6fBqzJLyJ\npKfnXQTBR2Li9fT09FBVVcXKlSvR6cR6mAGrm4E3G9DkRWNYknLS+4MBf0R2Vz7wNeZcdc0ZX3uw\n1832F2poruwjIcvMhkcKiEszXbDPNhbe2loG3tjMwJYtBDo6kOv1mNauJWrD9egXLZoWlU6CoSBN\ng01U91VT3VdNVV8VNX01o5amSTelkxedx7qsdaLcYvLIMGWguAD1TiUkxh1BEEX0aTmdTlwuKwQ8\npz6nNgp0FtDHgjEBEgrDz8OPodf04a0uBpTjOwdUEt6n6OrajE6XgdlUwpbtW1AoFJHoThAEbK+e\nALmMmFvyT+peDAWDvPXbX4Zl99Uzll0oJHB4awtlm+uRyWSs+FwexavSkF+k7svg4CADW7Yw8PIm\nPJWVoFBgWLaUhG99C9MVa5DrLt1FMj0BD7W22ojUqvuqOW47jico/iGr5Cpyo3NZnb6amZaZzLTM\npMBSIHVHSkxuQqFwpNULzt6wrHrBGd66rOHjQ8esEPSe4mQyUUZDYopKg+Q5oI8Ji2qEuHQj5KWY\n/DqZ/Hc4jvj9Nvpsu8nK/DIul4vDhw8zZ84cDAbxx851sBvviX6ib8hB+akldYRQiHf/+DjHubnD\nMQAAIABJREFU937C6nseYM5V68/omv3dLrY+U0Vn/QBZJXGsvD3/omRfCoKAe/9++l9+mcF33kXw\netHMnEniv30f8/r1l2RJL5ffRVVfFcd6j1HZV0m1tZqGwYZIl6RJZaLAUsAt+bdE5DYjaoY0l01i\n4gn4RkhrpLBOcczdB+Hv9UmoTWCIBX0cmFMhaY4oK0NcWFyxIyKxWDFSu0R7LiThjaC3dxsQIj5+\nLeXl5QQCgUhmZtDpZ2BLPeoME4ZFoyuFCILA1r/+kcodH7HstruZf+3Gz7yWIAgc+7iNTzadQKGU\nc9UXZ5G3MPGCj/0EensZeO01+l96GV9TE3KjkaibbiT65lsuqeQTb9BLTV8Nx6zHONZ7jGPWY9QP\n1EfklqBPoNBSyJWZV0bklmpMvWQ+v8Qkx+8GR/dwlHVSxPUpgXkHT3EiWVhMYVnF54v7hrjhY0Ny\nGxLaOHcbTmYk4Y2gt3cranUCOt1Mysp+S05ODgkJ4vy6gbcbCHmCxNyUd1JXZtlrL3H4/bdZuOFm\nFt30uc+8jsPm5aNnq2iu7CO9MIY19xRijLlwUZ0gCLj2lmF7/nnsH30EgQC6BfNJfugrmNetm/Jd\nlv6QnxO2ExyzHqOit4JKayW1tloCgrhytkVroTiumKsyr6IotoiiuCLidJdeBCsxgYRC4OkHZ4/4\nGJKZszv8PHzcGT5+qsQNuWpYVoZYiM4cHXlFXosb7jq8RKOv8UASXphQyIu1bweJiddRU3Mch8PB\nDTfcAICvzYGrvAvj8lRUSaPHcqp2bmPnP//OzGWrWHHHFz4zYmg82ssHz1QS9IdYdUc+RSsvXJQR\ncjoZ2LwZ2/PP4609gSI6Gss99xB9yy1oZow9V3CyIwgCbY42jvQc4XDPYSp6K6juq45M3DarzRTF\nFnFv8b0UxxZTFFdEov7CR8oS04CAb1hgIx9jyczVC6HAyeeQycOiShAllbpATOAwxA0fG5KbPg40\nJpC+q+OGJLwwNlsZwaCD+LgreOutA0RFRZGTkyOucbelHrleiXlNxqj3tFZW8O4fHydtVjHrHnr0\ntJPKg8EQe1+r5+D7zcSmGVn3QBExSRcmEcLb0IDtH/9g4JVXCTkcaGfNIvlnP8O8/hrk2vGtxnK+\nuPwujlmPRQR3pOdIpECyTqmjKLaIOwvvFCO32CLSTGmS3CROjd8Nji5RWkPbk0QWjsQ8A2OfQ6kD\nYzwY4sUxsORScd+YIG4jMosXuxulCGzSIgkvTK91K3K5FpmskPr6j1m9ejVyuRx3RS++hgGiN+Yg\n1w3/c/V3dfL6L39CVEISN3zrB6ctF2bv8/DekxV01g9StDKV5bfmnveq42K35V6sT/0V544doFJh\nXreOmLvuFJfcmQISEASBFnsLh3sOR+R23HacoBAEINOcybLUZZTElTAnYQ650bko5dJXdtoTCorj\nXI6u0TKzd50sN+8pJKazhGUVD0nF4f2hSCx+dFSmNkhR2CWC9OsRxmrdjiVmKYcPVwEwd+5chECI\n/rcbUCbqMSwcTlTxez288aufIiBw43d/eNpC0C3Vfbz7lwpCQYG1DxSRtyDxvO5TCASwv/ce1ief\nwlNZiSI2lrivP0zMbbehHLGKw2TEH/JTba3mQPcByrvKOdR9KDLXTa/UMzt+NvfPvp858XOYHTeb\nGO35rf4uMcXw2kfIKiwse+fJx5w9EP6folGoTaKoTEmixIyJ4nNjIhiTxCjNmCh2OUqZuNMSSXiA\n292G291Maurn+eCDg+Tm5hIVFYVjbwdBq4fYe4uQKcT/wxMEgfef+D96mhu56bs/JDpp7LXdBEHg\nyIetfLLpBDFJeq75ymyiT7GawpkQcrnof+VV+p5+Gn9bG+qsLJIe+xFRN9yAXDM5s7DcATdHe45S\n3l1OeVc5R3qORFYDSDelsyJtBaUJpcyJn0NOVI40iftSRBDE5I7BDrB3iAIb2n46GhuxBFIEuVKM\nsowJYEqGlNKwwEbKLLyvluZKSpweSXiAzbYbAI87E7u9mXXr1iEEQtg/bEGdYUJbMBxpHHxnC1U7\nt7HstrvJnrtgzPMF/SG2/aOG6l0dZM+J48r7ZqHWnts/ddDhwPbcc/Q98zeC/f3oSktJ/P73MK5Z\nM+lWJhjwDnCw+yAHug5Q3l1OpbWSQCiADBn5MflszN3IvMR5zE+YT7x+ckejEmeA1zFaYPb2Tz0P\nb8eqzqGNCkddCZA6X4zKIgIbEZXpYmCSfc8lpi6S8BCFp1JZqK11oVKpyM/Px7m/k+CAl5ib8yLj\nYV0NdWx/9ilmzL+MRTfeNua5PA4/b/7hCJ31Ayy4NovLrs0+p4LPQbsd23PPYX3mb4QGBjCuWkXs\nlx9EP2/86s59Fi6/i/Kucso6y9jbsZfqvmoEBJRyJcWxxXxh1heYlziP0oRSzGrzRN+uxJni94Cj\nc2x5DQ5JrRN89pPfqzKAOVmMxtIWiiIzJYe3KeFtEqim9tQYianJtBeeIAjYbLuJjl7Mnt1V5Ofn\no5IrsX4kRneavGhAHLd767e/QG82c/UpMjIHe91s/t1h7FYP675UTO78s1/RPGi30/fss/Q98zdC\ng4MYL7+cuK9+Fd3s4vP+rOeLN+jlSM8R9nTsoayjjIreCgJCAJVcRWlCKQ+VPsSCxAXMjpuNVjm1\nskOnDV6HKK3B1vC2HQZG7Ns7xKodn0ahHpZXYhHkXvkpmYW3Wul/bCQmL9NeeC5XA15fF6FgLi6X\njeLiYlwHuwkO+Ii5OT8S3W37+5P0dbRxy7//GJ3p5D/qnhY7W353mGAgxIZHS0nJjT6r+wh5vdie\ne57eJ54QI7o1a0TRFRddkM95LgRDQSqtlezp2MPezr0c6j6EN+hFLpNTHFvMvcX3clnSZZQmlKJT\nSv/HPuF47ScLbJTY2sbOWjTEgzkFYjIhY/GwvIYiNVOy2LUoZSpKTHGmvfAGBg8A0NyiRa1Wk5ub\ni/X/jqBKNkSiu/qD+zjywTssuP4mMmeXnnSOtuM23vzDETQ6JTc8Oh9LypkPngvBIAObN9Pzm98S\n6OjAsGIF8Y8+gq5oYkTX6+5lV/sudrbtZHf7bvq9/QDkx+Rza/6tLEpexPzE+ZjUF3cVB4lP4XWE\nRXaKyGywbexyVIaEsMyyIWu5uG9ODT9SRJmppGhcYnow7YU3OHgYhcJIdVU/+fn5BBscBLpcxNwq\nRnc+t4sP/vIHYtMyWH775096f0t1H2/9/gimOB0bvjHnjEuECYKAc+dOun/5K7w1NWiLikj5+c8w\njFhVfTzwh/wc6j7ErvZdfNL2CVV94rQMi9bCyrSVLE1ZyuLkxcTqYsf1vqYVgiBOgB5ohv4WUWQD\nLeH9ZvG523by+wwJEJUKsTmQvXJYZlEjZCbVUZSQiCAJb+AwGk0BLpeb/Px87DtakZvU6OeIWYQ7\n//ks9r5e7njsf1AoR8/daT5m5a0/HSU6QccNj85FZ1Kf0TW9DQ10/eznOHfsQJWeTur//grT1VeP\nW9Zlr7uXj1s/ZnvLdvZ27sXpd6KUKZmTMIdH5j3CspRlFFgKkMuk7LgLQtAfjsjCMhsSWX/L8LFP\nZzKqjRCVDtHpYvJHVHr4EY7OTMmgPLPvm4SEhMi0Fl4w6MbhrEYQrkYmk5EZk4a99hjmtZnIlHLa\nj1dz8N0tlK5dT0p+4aj3Nh2z8vYfjxKTrGfDI6XojJ/94xNyOun905+xPvMMco2GhO99F8uddyJT\nX9wfLkEQqB+o56OWj9jWso0jPUcQEEgyJHFN9jUsT1nOZcmXSd2U50rQL0rL1gj9TWBrGhGhtYrp\n+p9eusUQLwossQjyrx6WW1S6uP6YNGYmIXHBmdbCG7RXIAhB2ts1pKWlwbEBkINhQRKhUJAPnvw9\nxhgLy2//wqj3tZ/o5+0/ibK74dG5aA2nr9ogCAL2d96h67/+m0BXF1EbN5LwrW9e1MoogVCAg90H\nI5JrsbcAMCt2Fg+VPsSa9DXkx+RPiRJkE44giNU9bE1hqTWO2G8Sk0FGVv6QK8UuxagMyF4xLLHo\ndPFYVKqUli8hMQFMb+ENHgKguUnOiuW5OHd2oS2woDCrObL1XXqaGrj2ke+g0Q9XSOltdfDm749g\nsmjZ8I3Sz5Sdv7OTzv/8EY5t29DOmkXq479GP3fuRfk8/pCfso4y3mt6jw+bP6Tf249KrmJR8iLu\nLbqXlWkrSTIkXZRrT3m8juHorD8ss5FS87tGtzcmiku5pC+GkkxxPyZLzHQ0p0oFhCUkJiHTXHhH\nkcvj8ft1pMnjCdmtGBYm4XO7+OTFZ0nOn0nBkhWR9gM9bjb/9hAqjULsxjzNmJ0QCtH/r5fo/sUv\nEEIhEr//PWLuvhuZ4sL+EA5J7t3Gd/mw5UMGvAPolXpWp6/myswrWZqyFINKKrkEgLsf+urHfjh7\nRrdVG0WBWWZAzhpRZENSi84A9bmXiZOQkJgYprXwnM7j+H2JqNVqjHVBAiY12gILn7z0HK6BfjZ+\n+z8iXX5ed4A3fy/Os7vp/83HZDl1NqavtZWO7/8brn370C9ZTPJjj6FOT79g9x0MBdnbuZd3Gt6J\nSM6gMrA6fTVrM9eyLHUZGsU0zc5z20SBWUcKrU7cuqyj25pTRaEVXCOm7ccMCS1LXOZF6u6VkLik\nmLbCC4W8uFz12GzzSU9Nx3d8AOOSFFz2fsrffI2Zy1aRnFcQbivw3pPHGOh2c/0jpaedZzeweTOd\nP3oMgOSf/oSom266YONkNX01bK7bzFsNb9Hj7pm+kvM6wFoLvSdEmVnrhsX26fR9cxrEzoDC68GS\nIwrOMgMs2dI4moTENGPChCeTydKBvwNJQAh4QhCE34zX9Z2uBgQhSG+vmsLYOAgK6OfEs+uNFwn6\n/Sy99c5I292vnKD5mJVVdxaQVjD2kjVBu53OHz3G4JYt6ObPJ+W//xt1Wup532eXs4u3Gt5ic/1m\nam21KGVKlqct5/oZ17MqfdWlKzlBECdT99aGH8fDkqsVj0eQiQkhlhkwa6M4J21IajFZktQkJCQi\nTGSEFwC+JQjCAZlMZgLKZTLZ+4IgVI7HxZ2OGnHrjCFeoUcRrcFn9HP4/bcpXL6KmGRRVnUHujn0\nQQuzV6dRvHJsgbmPVtD26KP4OzuJf+QbxD744HmN1flDfra1bOPl4y+zu303AgIlcSX826J/4+qs\nqy+tdeL8bjFC6z0uyswallvvidHLxWjMEJsLWSsgLg/i8sVtTLZUKURCQuKMmDDhCYLQAXSE9+0y\nmawKSAXGRXgORw2CoMDvt2BukaFbFk/5m68R9PtZdNPtgFgM+sNnq0nIMrPsltwxz9O/aROdP3oM\nRVwsWc8/h6705NJjZ0rLYAubajfx2onXsHqsJOoT+VLJl7h+xvVkRWWd83knBX63KLLuauipgu7w\no78ZEIbbRWWIIpu3FOJyw2LLF7MipTE1CYkphSAIeAMhnN4ALl8w/BjeL0g0kRE7fglgk2IMTyaT\nZQFzgb1jvPYg8CBARkbGBbum03UCv99CojkRhV2GKs/IkZ++Q/7iZVhSUgkGQrz75DEA1j1QhEI5\nuupIyOej6yc/pf9f/8KwdAkpv/oVypizj7z8IT8fNn/Iy8dfZk/HHuQyOSvTVnJr/q0sS1k29RZF\nDfjAemK01LqrwNYwPPlarhKlljofSu8cjtgsOVL2o4TEBBEKCbj9QZy+AC5vEEdYUkPPxW0Apy8Y\nEZjzU23E94yWW0g49TV/tKGILyzNGrfPOOHCk8lkRmAT8KggCCdVvxUE4QngCYAFCxac5p/u7HC5\nGrHbdcQJZuR6Jcfr9uJzu5i3/gYAyrY00N04yNUPFmOOGz0OFLDZaP3q13AfPEjsgw8S/8g3zroL\nc8A7wEvHX+If1f+g29VNsiGZr5V+jRtzbyTRkHihPubFQxDEaiKdR6GzAroroadalF0oILaRKcSx\ntMQimH0LJBRCfKE4zqY4/fxFCQmJUxMKCbj8wVEC+rR8nL5hKTnD7Vy+AE7v8NY54rnLF/zsC4dR\nymUYNEoMagX6oa1aSUq0Fr1aiUGjQKcKb9UKDGplZKtXK8IPJakx4zvGPqHCk8lkKkTZPS8Iwivj\ndV1BCOF2t+B25ZLYr0GTG83Bd/9GUk4eyXkF9DTbOfheMzOXJpMzb/Sadr6mJpoffJBARyepj/8a\n89VXn9W16/vreb7qed6oewNP0MOi5EX8x+L/YEXqiskbzQW8osyG5NZ5FLqOgmdoqRmZmCCSUAgz\nrxWlljATYvOk8TUJCSAQDOH0BnH4Ajg8ARxhCTm8w/sRKXkDOIak5BuWmigpUVBu/3nIKbyfEq3G\noBHFY1ArxONDzyPHxX2DRhSVQa1Er1GgVsinZJWmiczSlAFPAVWCIPzveF7b6+tGEHx4PEZi3Xoc\nukFs7a2sf/hbhEICW/9ehc6kYtnNo8ft3IcO0fLQV0EQyHjmGfTzzrxiyoGuA/zl6F/Y2bYTtVzN\ndTnXcVfhXeTH5F/oj3d+uPqg88houfXWDEdtKr0YsRXfDInFkFQCibNALU1ul7i08AaCYhTkDWD3\nBHD6woLynCwrR7g7z+kNjLn1+EOffUGG5WQMC2YoeorWqyPiGYqmhqQ01Nbwqa14XIlaKRWBH2Ii\nI7xlwOeBozKZ7FD42L8JgvDWxb6w2y3WlfT7ookS9Byt+xidOYr8Jcs5+F4z1lYH13xl9qiyYc49\ne2h56Kso4+PJeOLPqLOyPvM6giCwu303Txx9gvKucixaC18r/Rq3FdyGRWu5WB/vzPEMQPshaD8Y\nfhwIJ5GEMSVD0mwouHpYbpZsqWyWxKTFGwji8IiCcnjHjqQc3uCwtHxjC8zpDeILnpmkdCpFWFKK\niKySzFpxXys+N6iH9ofbGMPR1tC+XqNAo5T+ti4mE5mluROYkJjY7W4CQOdLQZmo4/j+T5hz1Xrc\n9iDlbzeSMzeeGaXDhZ0dOz+h9WtfQ52RTsbTT6OMizvt+QVBYFvLNv585M8csx4jQZ/Adxd+l5vz\nb564lcG9DjFyi8jtoDjeNkR0JqTMgwX3Q/IcUXSG039OCYkLRSgk4PSJohJl5WfQE4jIy+7x4/AO\nv273+CNSG9q3ewP4AmcmKaNmuKvOFBZPukEf2R8S2CgpaZUnycqgVqBUSBHUVGHCk1YmAre7GUGQ\nobfH4YgbJBgIMGvVFex+tQ4hBEtHdGU6tm+n9eGvo87JIeOvT6G0nD4y29+5n8cPPM7hnsOkGdP4\n4ZIfsiFnA2rFOK5dFgqJ89la9kJLGbTuF7slh7IkzamQMhfm3CFuU+aKpbQkJM4BbyAoymekjLwj\nRBWW0ViiikRjvgDCZ6SkyWWiqExaFSatEpNWSZxRTXacAZNWFJI5/NqQkEYKzBQWll6lQC6feuNP\nEufPtBSew9GIx2Mgym+isf0wcemZhEKxHC87wIL1WZGsTNf+/bR+4xE0eXlk/PUpFNHRpzxndV81\nvznwG3a27SRBl8APl/yQG3JvQCUfh2xEr12UWus+UXKt+4YTSnQx4gKiRRtFsSWXgmkKZIFKjAtD\n86QG3H4G3X4GPX4G3YHwVoyyho6LbYZfG4q2zqTrT6uSY9SoMIdFZdQqiTcaI/smrQpTWEomrSp8\nTIlZq8SoESWmVyumZKKExORhmgqvCa/XQJKg53DjPubfeiNlmxvQmdXMXSvO9fPU1NDy0FdRpaSQ\n/tSTp5Rdr7uXx8sf5/W61zGrzXxz/je5Y+YdaJUXMTtxsAOaPhEfLWXilAAhBMjETMlZGyF9EaRf\nJlYnkX4kLmk8/uBZiWrQE8A+Qm6fJSytSo5Zq8KsE4VlMajJjDWE5TUcbY2Mvoya4WjLoJESJyQm\nB9NSeD5fNz6vAZNahztox5xYxP53u1h2Sy5qrRJ/ezstD3wJuV5PxpN/GXNCuT/k54WqF/jT4T/h\nCXq4r+g+Hih5ALPafOFveKAVGj+Bpp3QuFMslAxiua20BTDzOlFuaQtAG3Xhry9x0fH4gwy6/fS7\n/fS7/PS7fAy4RVn1u8LbUVGYKK4Bt/8zx63UCrkoK50ooSidigyLHrNWGZbY8GtDUjPrxHYmrVJK\npJC4ZJh2whMEgWCwD78vjlDAQXRSMjV7vehMKopWphJyu2l9+OuE3G6y/vECqtST62eWd5Xz2O7H\nqB+oZ3nqcr678LsXtvTXYAfUbxPl1rRTXIQURJllLhMTS7KWi4klUsbkpEEQBOzeAAOukZLyRfZF\neYnPh+Ql7vtOm7Yul0G0Xo1ZqyRKJ0opJVp3SlENSW3oNa1K+o5ISMA0FF4gMIhMFkDujaLLWkdy\n0VwaKvpZenMuSpWc9n/7//BUVZH2h9+jycsb9V6X38XjBx7nH9X/INWYyu/W/I5VaavOf1zB54Km\nXVD3IdR/JHZRgjj+lrkMFn1F3CYWSYIbBwRBwOUL0uf0YXP5sLn82EbsD7h8kUhspMgG3P7TllHS\nquRE69RE64ejrJI0FdF6NVHhiCparxrVJkqvwqhWSkkWEhIXgGknPK+vGwClN5o+dwcadxEqrYKi\nFSnYXniBwS1biH/0EUyXXz7qfXs69vCfu/6Tdkc7dxfezdfnfh296hzrPgoCdFXAia2i5Jp3Q9AH\nCg1kLhWzJ3Muh4QikEtjH+eDIAg4fUFsTt8IgfmwOf3YXOKxfpf/pNdONa4lk4FZqxoWkk5FukVP\n9AhhiVv1CIGJEZgUaUlITCzTTng+ryg8tTeGXlqwNeiYvSoZobWR7v/5BYaVK4j98pcj7f0hP787\n+DuerniaLHMWf7vmb8xNOPMKKxECXmjcATVvQ807MNgqHk8ogssehJw1ouyk9dtOiy8Qos/po9fh\npc/pw+r00ucUI7A+l4/+sMSGhGZz+fAHxw67hroKY/QqYvRq0i165qRFE21QYdGridGriTGosRhE\ngVn0asw6FQop2pKQmJJMO+E5ne0AqL1RKKLd+D1yipYl0vb1+5AbDKT87GeRLspWeyvf/fi7HOk9\nwm35t/Hthd8+u+xLpxVq34Oat8RIzucQS3PlrIHV34PcK8GcfDE+5pQhGBKwuXxYHaK8rA4fVocX\nq9MnPhzhY+H9QU9gzPPIZUQEFaNXkRmrpzQ9+iRhDb1uMagxa1VSV6GExDRi2glv0C5GViG7Cr8v\nibSZMQQ2PYO3qoq0P/whUkVlZ9tOvrP9OwD8atWvWJu19swu4OqDqjfg2KvQ8LE4XcCYBLNvhYJr\nIHvlJR/FuX1Beh1euu0eeuxeeh2i0PqcXnpHSKwvHJWNNeF4SGCxRjWxBg1FKWZiDWpijZrIsTij\nGotB3DdppXEuCQmJ0zPthOdydhMKyRHcAfy+JGZkK+j7wVOYN1yPac3lCILAc1XP8cv9vyQ/Jp9f\nr/41aaa005/UbYPqN6HiFTG7UgiKK3EvexQKrxcne0/xsbhgSKDP6YtIrMfupcchbrvDz3vDW7t3\n7CjMrFUSFxZWTryRhdlq4kZIzGJQi68b1ETr1VLXoYSExAVl2gnP7eolEFAT9LpQapLQv/57Alot\nid/5DoFQgJ/s+QmbajdxRcYV/Gz5z06dmBIMQN1WOPicOC4X8kN0Biz9OhTdKNajnAITvr2BIN2D\nXjoHPXQNDstsSGJDYrM6vGNmIBo1SuJNGuJNGgpTzKw0aiLP400a4sPPY/RqafKxhITEhDLthOfz\n2vD7NfgDQdJSZfie/5jEH/yAUIyZ72z/f2xt3sqXZn+Jh+c+jFw2xg90dzUceh6OvAiOLtDHwsIH\nxC7L1HmTRnKCINDv8tM56BFlNuCJSK1zwEPnoJeuQQ99Tt9J71XKZRFhJUdpKUmLijxPiIhMS5xJ\njV497b5CEhISU5Rp92sVDAwQCmhxBZWkHX0bdXY2mpuv5+GtD7O7Yzffu+x73FV41+g3BXziuFzZ\nX6BlD8iVkLcO5t4FuVeBchwLQyNWlu91eGnrd9Pe76FjwC2KbNA7SmzeMSpwxBnVJJq1pERpmZsR\nTZJZS5JZS2KUlkSzhgSTlmidlMwhISFx6THthCcIDgS/Dhd6TJWvEvNfP+TRHd+irLOMx5Y+xo15\nNw43HuyA8meg/GkxmrPMgLU/hZLPgTH+lNc4X1y+AO39Htr73ZFHW7+Htn5XRHCfTrXXquSiuMzD\nIks0a0mKErdDMpO6FSUkJKYr0094Mg8EotEHwTAjlZ/oP2JPyx5+vOzHbMzdKDbqPQGf/BoOvyiu\n9J23dniu3AVIPrF7/LT0uWnuc9Fqc9FqG5KauLW5/KPay2WQZNaSEq2jND2aa0uSSYnWkRqtJTlK\nR0qUDrNOKVWSl5CQkDgN0054crkPIaAmrquRXVfG8X7LB3x7wbdF2XUchh3/C5Wvg1IDC+6DxQ+J\nkd1ZEAiG6Bjw0NLnonnEo8XmpqXPddK4mVGjJDVaR0q0ltL06LDMdKTG6EiJ1pFo0kiLTEpISEic\nJ9NPeEo/+FVEBzr4oaWKL8y6l3sSFsM/74LqLeIKBCu+CYseOm23pT8YotXmpr7HQUOvk/peJ81W\nUWzt/W4CI1IalXIZqTE6Mix6ri5OIsOiJz1GL24tOqL14zsGKCEhIfH/t3e3MVacZRjH/xewdGV3\n2SJb6FapLOnSFEtqkQA2RmnqS0sifGkM1UZriMRW+4E2JkY/1LR+sjEmTZpUjMSX+IIaoxulaaLW\nYrA0JSHFtkqLW8RNawFRsCBvy+2HmWzawmGnye48e+a5fglhZs+w5+bOOXvtPPOcZ3KUVeCNjp5m\n2rRR4mwH2wdeYPkV72XTy3+HbSuhowtWfxlWfW7sFjsRwavHTjF8uAi1lw4dL/4+fJwDR068IdTm\nzOrgyrldXLfgUj52XX8ZZkWoXT6702doZmaJZRV4J8u7gMfpDva8Zzqbd/+W6WdPESs2MrL08/z1\n2CW88OQhXnx1mBcPvsbwoeP878zo2L/v7JjGQF831/TPZs3Sfgb6uhi4rIuBuV3M6fIPEV0+AAAG\nA0lEQVRZmpnZVJZV4B15pVhWLE5PY9OZ/YzMXMY9szaw40/dnHziz2PHXdHbyVXze1gx8HYWXdbN\nor4uBvq6uHx2p6frm5m1qawC7+nfbKV3KZwbPceWY3exp+v9DM7p4ZNX9bB4fjeD83sYnNdNT2dH\n6lLNzGyCZRV47161lmd2/JPeBTfy4N2foPdtDjYzs1xkFXhXr7yBq1fekLoMMzNLwFMHzcwsCw48\nMzPLggPPzMyy4MAzM7MsOPDMzCwLDjwzM8uCA8/MzLLgwDMzsyw48MzMLAsOPDMzy4IDz8zMsuDA\nMzOzLDjwzMwsCw48MzPLggPPzMyy4MAzM7MsOPDMzCwLDjwzM8uCA8/MzLLgwDMzsyw48MzMLAtJ\nA0/SzZL2Ston6UspazEzs2ZLFniSpgMPA7cAS4DbJC1JVY+ZmTVbyjO8FcC+iBiOiNPAT4B1Cesx\nM7MGSxl47wD+8br9kfJrZmZmE25GwufWBb4W5x0kbQQ2lruvSdo7Ac/dBxyegO/TRO5Na+5Na+5N\na+5NaxPRm3dVPTBl4I0AC163/07g5TcfFBGbgc0T+cSSdkXE8on8nk3h3rTm3rTm3rTm3rRWd29S\nDmk+DQxKGpA0E1gPDCWsx8zMGizZGV5EnJX0BeAxYDqwJSKeS1WPmZk1W8ohTSJiG7AtwVNP6BBp\nw7g3rbk3rbk3rbk3rdXaG0WcN0/EzMyscby0mJmZZaGxgTfesmWSLpG0tXz8KUkL668yjQq9uUfS\n85L2SPqdpMrTfttd1eXuJN0qKSRlM/uuSm8kfbx87Twn6Ud115hKhffUlZIel7S7fF+tSVFnCpK2\nSDoo6dkWj0vSQ2Xv9khaNmnFRETj/lBMgvkbsAiYCTwDLHnTMXcBj5Tb64GtqeueQr25EZhVbt/p\n3px3XA+wHdgJLE9d91TpDTAI7AbmlPvzUtc9hXqzGbiz3F4C7E9dd439+QCwDHi2xeNrgEcpPpu9\nCnhqsmpp6hlelWXL1gHfK7d/Dtwk6UIfhm+acXsTEY9HxIlydyfFZyRzUHW5uweArwMn6ywusSq9\n+SzwcET8GyAiDtZcYypVehPA7HK7lwt85ripImI7cOQih6wDvh+FncClkvono5amBl6VZcvGjomI\ns8BRYG4t1aX1Vpd020Dx21cOxu2NpOuBBRHx6zoLmwKqvG4WA4sl7ZC0U9LNtVWXVpXefBW4XdII\nxcz0u+sprS3Utsxk0o8lTKIqy5ZVWtqsgSr/vyXdDiwHPjipFU0dF+2NpGnAN4E76ipoCqnyuplB\nMay5mmJU4I+Sro2I/0xybalV6c1twHcj4huS3gf8oOzNuckvb8qr7WdxU8/wqixbNnaMpBkUwwwX\nO+1uikpLukn6EPAVYG1EnKqpttTG600PcC3wB0n7Ka43DGUycaXqe+pXEXEmIl4C9lIEYNNV6c0G\n4KcAEfEk0EmxjqRV/Jk0EZoaeFWWLRsCPl1u3wr8PsorqA03bm/KYbtvUYRdLtdhYJzeRMTRiOiL\niIURsZDi+ubaiNiVptxaVXlP/ZJiwhOS+iiGOIdrrTKNKr05ANwEIOkaisA7VGuVU9cQ8KlytuYq\n4GhEvDIZT9TIIc1osWyZpPuBXRExBHyHYlhhH8WZ3fp0FdenYm8eBLqBn5XzeA5ExNpkRdekYm+y\nVLE3jwEfkfQ8MAp8MSL+la7qelTszb3AtyVtohiuuyOTX7CR9GOKYe6+8hrmfUAHQEQ8QnFNcw2w\nDzgBfGbSasmk52ZmlrmmDmmamZm9gQPPzMyy4MAzM7MsOPDMzCwLDjwzM8uCA8/MzLLgwDMzsyw4\n8MzaSHlPtQ+X21+T9FDqmszaRSNXWjFrsPuA+yXNA64HGr8CjtlE8UorZm1G0hMUS7+tjoj/pq7H\nrF14SNOsjUhaCvQDpxx2Zm+NA8+sTZR3gf4hxR2ij0v6aOKSzNqKA8+sDUiaBfwCuDci/gI8QHEX\nbTOryNfwzMwsCz7DMzOzLDjwzMwsCw48MzPLggPPzMyy4MAzM7MsOPDMzCwLDjwzM8uCA8/MzLLw\nf5x9ZmdQ28pjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10fc28750>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot example solutions\n",
    "for Ts in Tsol[::10]:\n",
    "    plt.plot(np.linspace(0,1,N), Ts)\n",
    "plt.xlabel(r\"$x$\")\n",
    "plt.ylabel(r\"$T$\")\n",
    "plt.ylim([0,10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No re-use of solution:\n",
      "[ 9 10 10 11 11 11 12 12 13 13 13 14 14 15 15 16 16 16 16 17 18 18 18 18 19\n",
      " 19 20 20 21 21 21 21 21 23 23 23 24 24 24 24 24 26 26 26 27 27 27 27 29 29\n",
      " 29 29 29 30 30 31 32 32 32 32 32 33 34 34 35 35 35 35 35 37 37 38 38 38 38\n",
      " 38 38 40 41 41 41]\n",
      "With re-use of solution:\n",
      "[ 9  9  9 10 10 10 10 10 11 11 11 12 12 12 12 13 13 13 13 13 13 13 14 15 15\n",
      " 15 15 15 16 16 16 16 16 16 16 16 18 18 18 18 18 19 19 19 19 19 19 19 19 21\n",
      " 21 21 21 21 22 22 22 22 22 22 22 22 24 24 24 24 24 25 25 25 25 25 25 25 25\n",
      " 27 27 27 27 27 28]\n",
      "Total iterations with / without re-use: 1472 / 2031\n"
     ]
    }
   ],
   "source": [
    "# Same as previous approach, but now take advantage of the previous result\n",
    "# as the starting guess for the next value\n",
    "Tsol = np.zeros((Tf.size,N))\n",
    "it_save = np.zeros(Tf.size, dtype=int)\n",
    "for i in range(Tf.size):\n",
    "    if i == 0:\n",
    "        T0 = None\n",
    "    else:\n",
    "        T0 = T\n",
    "    T, itcount = Tsolve(Tf[i], N, K, Kp, T0=T0)\n",
    "    Tsol[i,:] = T\n",
    "    it_save[i] = itcount\n",
    "    \n",
    "# Compare iteration counts\n",
    "print(\"No re-use of solution:\")\n",
    "print(it)\n",
    "print(\"With re-use of solution:\")\n",
    "print(it_save)\n",
    "print(\"Total iterations with / without re-use: {:d} / {:d}\"\n",
    "      .format(np.sum(it_save), np.sum(it)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus we have solved the problem with relaxation, and we see that re-using saved solutions reduces the workload by about 25% compared to starting from scratch."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Further Improvements\n",
    "\n",
    "A further improvement that could be made to our scheme, but that we shall not implement here, is adaptive gridding. Instead of just keeping $N$ fixed, we could try a guessed value of $N$, try again with $2N$, and compare the results. If they agree within some tolerance, we're done. If not, we increase $N$ further until we get a solution that we're happy with. This can be made even more sophisticated by making the grid non-uniform in space, for example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}