notebooks/lao-3-3-solved.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Loops & Orbits &mdash; Week 3 &mdash; Day 3 &mdash; Jupyter Notebook\n",
    "\n",
    "## Ellipses from Newton's Cannon with Second-Order Runge-Kutta\n",
    "\n",
    "Pronounce the German names well: \"roon-geh koot-tuh\" &mdash; not \"runj kutah\" which is admittedly how it looks.\n",
    "\n",
    "**Execute the imports.** There is nothing here for you to change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Important math functions and the constant pi:\n",
    "from math import sin, cos, atan2, sqrt, pi\n",
    "\n",
    "# We will start using the high-performance numpy library for dealing with arrays of floats:\n",
    "import numpy as np\n",
    "\n",
    "# The following import statement makes the plotting library available to us. There is also a mysterious\n",
    "# statement beginning with % to work around a known Jupyter bug: https://github.com/jupyter/notebook/issues/3670\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Use unit testing framework to test each function we write:\n",
    "import unittest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, define all the constants. *Our numbers in meters were getting really huge. All of the constants are now in kilometers instead of meters.* **Execute these assignments.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The following two parameters are get us the strength of gravity at various heights:\n",
    "g = 0.00981                    # 0.00981 km/s**2 at Earth surface\n",
    "radius_of_earth = 6371.0 # Earth radius in km\n",
    "\n",
    "# The following two variables define the cannonball launch conditions:\n",
    "initial_position = [0.0, radius_of_earth + 500.0]  # cannon is on a 500 km high mountain\n",
    "initial_velocity = [8.982, 0.0]                     # a good range of initial x-velocities to try is 7.6 to 9.2 km/s\n",
    "# The value of 8.982 km/s was tweaked to make exactly one elliptical orbit in 200 minutes.\n",
    "\n",
    "# The following two parameters establish the simulation time step and its maximum total duration:\n",
    "delta_t = 120.0              # simulation time step in seconds\n",
    "simulation_time_steps = 100  # 200 minutes worth for a delta_t of 120.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions for Working with Vectors\n",
    "\n",
    "All of these functions you implemented previously. **Just execute this cell to get our vector functions.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "test_angle_from_vector (__main__.VectorTests) ... ok\n",
      "test_length_from_vector (__main__.VectorTests) ... ok\n",
      "test_vector_from_length_and_angle (__main__.VectorTests) ... ok\n",
      "\n",
      "----------------------------------------------------------------------\n",
      "Ran 3 tests in 0.003s\n",
      "\n",
      "OK\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<unittest.runner.TextTestResult run=3 errors=0 failures=0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# computes horizontal and vertical components of a vector and returns them as a tuple\n",
    "def vector_from_length_and_angle(length: float, angle: float) -> np.ndarray:\n",
    "    # we are working in degrees -- python's are expecting radians -- convert degrees to radian\n",
    "    angle_in_radians = angle * pi / 180.0\n",
    "    x_component = length * cos(angle_in_radians)\n",
    "    y_component = length * sin(angle_in_radians)\n",
    "    return np.array([x_component, y_component])\n",
    "\n",
    "# get angle from components using atan2 version of arctangent function -- result is converted to degrees\n",
    "def angle_from_vector(vector: np.ndarray) -> float:\n",
    "    # use the arctangent function\n",
    "    angle_in_radians = atan2(vector[1], vector[0])  \n",
    "    # we are working in degrees -- python's functions return radians -- convert radians to degrees\n",
    "    angle = angle_in_radians * 180.0 / pi\n",
    "    # return the result\n",
    "    return angle\n",
    "\n",
    "# get length from components using Pythagorean theorem\n",
    "def length_from_vector(vector: np.ndarray) -> float:\n",
    "    length_squared = np.sum(vector**2)\n",
    "    return sqrt(length_squared)\n",
    "\n",
    "### Unit tests ###\n",
    "\n",
    "class VectorTests(unittest.TestCase):\n",
    "    \n",
    "    def test_vector_from_length_and_angle(self):\n",
    "        vector = vector_from_length_and_angle(1000.0, 30.0)\n",
    "        expected_vector = np.array([866.02540378, 500.0])\n",
    "        self.assertTrue(np.allclose(vector, expected_vector),\n",
    "                        msg=\"Expected {}, but got {}\".format(expected_vector, vector))\n",
    "        \n",
    "    def test_angle_from_vector(self):\n",
    "        angle = angle_from_vector(np.array([866.02540378, 500.0]))\n",
    "        self.assertAlmostEqual(angle, 30.0)\n",
    "        \n",
    "    def test_length_from_vector(self):\n",
    "        length = length_from_vector(np.array([5.0, 12.0]))\n",
    "        self.assertAlmostEqual(length, 13.0)\n",
    "\n",
    "testSuite = unittest.TestLoader().loadTestsFromName(\"__main__.VectorTests\")\n",
    "testRunner = unittest.TextTestRunner(verbosity=2)\n",
    "testRunner.run(testSuite)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions to Get Acceleration of Gravity\n",
    "\n",
    "All of these functions for gravity you implemented in Notebook 3-1. **Just execute this cell to get the gravity functions.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "test_acceleration_of_gravity_4g_up (__main__.GravityTests) ... ok\n",
      "test_acceleration_of_gravity_g_left (__main__.GravityTests) ... ok\n",
      "test_direction_of_gravity_left (__main__.GravityTests) ... ok\n",
      "test_direction_of_gravity_up (__main__.GravityTests) ... ok\n",
      "test_strength_of_gravity (__main__.GravityTests) ... ok\n",
      "test_strength_of_gravity_high_up (__main__.GravityTests) ... ok\n",
      "\n",
      "----------------------------------------------------------------------\n",
      "Ran 6 tests in 0.005s\n",
      "\n",
      "OK\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<unittest.runner.TextTestResult run=6 errors=0 failures=0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def strength_of_gravity(position):\n",
    "    # this function encodes the strength of gravity as a function of distance from the center of the Earth\n",
    "    radius = length_from_vector(position)\n",
    "    strength = g * radius_of_earth**2 / radius**2\n",
    "    return strength\n",
    "\n",
    "def direction_of_gravity(position):\n",
    "    # this function encodes the direction of gravity (the angle)\n",
    "    # gravity is attractive -- it always points toward the center of the Earth\n",
    "    direction = angle_from_vector(position) + 180.0\n",
    "    return direction\n",
    "\n",
    "def acceleration_of_gravity(position):\n",
    "    # using the strength and direction functions you have just implemented compute and\n",
    "    # return a 2x1 array for the acceleration of gravity\n",
    "    strength = strength_of_gravity(position)\n",
    "    direction = direction_of_gravity(position)\n",
    "    acceleration = vector_from_length_and_angle(strength, direction)\n",
    "    return acceleration\n",
    "\n",
    "### Unit tests ###\n",
    "\n",
    "class GravityTests(unittest.TestCase):\n",
    "\n",
    "    # tests of strength\n",
    "\n",
    "    def test_strength_of_gravity(self):\n",
    "        position = np.array([radius_of_earth, 0.0])\n",
    "        strength = strength_of_gravity(position)\n",
    "        self.assertAlmostEqual(strength, g)\n",
    "\n",
    "    def test_strength_of_gravity_high_up(self):\n",
    "        position = np.array([0.0, 2.0 * radius_of_earth])\n",
    "        strength = strength_of_gravity(position)\n",
    "        self.assertAlmostEqual(strength, 0.25 * g)\n",
    "        \n",
    "    # tests of direction\n",
    "        \n",
    "    def test_direction_of_gravity_left(self):\n",
    "        position = np.array([100.0, 0.0])\n",
    "        direction = direction_of_gravity(position)\n",
    "        self.assertEqual(direction, 180.0)\n",
    "\n",
    "    def test_direction_of_gravity_up(self):\n",
    "        position = np.array([0.0, -100.0])\n",
    "        direction = direction_of_gravity(position)\n",
    "        self.assertEqual(direction, 90.0)\n",
    "        \n",
    "    # tests of acceleration\n",
    "\n",
    "    def test_acceleration_of_gravity_g_left(self):\n",
    "        position = np.array([radius_of_earth, 0.0])\n",
    "        acceleration = acceleration_of_gravity(position)\n",
    "        expected_acceleration = np.array([-g, 0.0])\n",
    "        self.assertTrue(np.allclose(acceleration, expected_acceleration),\n",
    "                        msg=\"Expected {}, but got {}\".format(expected_acceleration, acceleration))\n",
    "\n",
    "    def test_acceleration_of_gravity_4g_up(self):\n",
    "        position = np.array([0.0, -radius_of_earth / 2.0])\n",
    "        acceleration = acceleration_of_gravity(position)\n",
    "        expected_acceleration = np.array([0.0, 4.0 * g])\n",
    "        self.assertTrue(np.allclose(acceleration, expected_acceleration),\n",
    "                        msg=\"Expected {}, but got {}\".format(expected_acceleration, acceleration))\n",
    "\n",
    "testSuite = unittest.TestLoader().loadTestsFromName(\"__main__.GravityTests\")\n",
    "testRunner = unittest.TextTestRunner(verbosity=2)\n",
    "testRunner.run(testSuite)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The While Loop That Does the Work\n",
    "\n",
    "Here is where your work lies.\n",
    "\n",
    "You have learned about Kepler's laws for ellipses. Euler and Euler-Cromer have too much accumulated\n",
    "error to work well with ellipses.\n",
    "\n",
    "**Below you will be implementing second-order Runge-Kutta.** Look for the five TODO comments in the code. That's where your work is. I have tried to make second-order Runge-Kutta understandable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the x and y velocities\n",
    "velocities = np.zeros((simulation_time_steps, 2))\n",
    "velocities[0] = initial_velocity\n",
    "\n",
    "# Initialize the x and y positions\n",
    "positions = np.zeros((simulation_time_steps, 2))\n",
    "positions[0] = initial_position\n",
    "\n",
    "# Initialize the times\n",
    "times = np.zeros((simulation_time_steps,))\n",
    "\n",
    "# The 0th elements of velocities, positions, and times are already filled in, so start with 1\n",
    "for i in range(1, simulation_time_steps):\n",
    "    #\n",
    "    # get all the before values\n",
    "    #\n",
    "    # velocities\n",
    "    before_velocity = velocities[i - 1]\n",
    "    # positions\n",
    "    before_position = positions[i - 1]\n",
    "    # time\n",
    "    before_time = times[i - 1]\n",
    "    #\n",
    "    # fundamental change for 2nd-order Runge-Kutta -- first estimate mid_position!!\n",
    "    #\n",
    "    mid_position = before_position + before_velocity * delta_t / 2.0\n",
    "    #\n",
    "    # calculate mid_acceleration -- using the estimated mid_position!\n",
    "    # \n",
    "    mid_acceleration = acceleration_of_gravity(mid_position)\n",
    "    #\n",
    "    # calculate after_velocity from before_velocity and mid_acceleration\n",
    "    #\n",
    "    after_velocity = before_velocity + mid_acceleration * delta_t\n",
    "    #\n",
    "    # calculate after_position using the democratic combination of after_velocity and before_velocity\n",
    "    #\n",
    "    after_position = before_position + (before_velocity + after_velocity) / 2.0 * delta_t\n",
    "    #\n",
    "    # update time\n",
    "    #\n",
    "    after_time = before_time + delta_t\n",
    "    #\n",
    "    # assign the after values into their lists\n",
    "    #\n",
    "    velocities[i] = after_velocity\n",
    "    positions[i] = after_position\n",
    "    times[i] = after_time\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diagram\n",
    "\n",
    "**Execute this to make a graph that looks a little like [the diagram Newton made](https://upload.wikimedia.org/wikipedia/commons/thumb/2/2f/Newton%27s_Principia_%281846%29.djvu/page519-1024px-Newton%27s_Principia_%281846%29.djvu.jpg).** Our ellipses are better than Newton's!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5AAAANcCAYAAAApZkDrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzde3xcZYH/8e/TJG0naelQqAsZkBbFcIsQGS5LxeWihiKXGIUVL4sLyMqPRX8qESLL6nohlbC4eOOnIgjsijdCAEEjUBQUFVIjhAKBcmfCcnVom5zcn98fTwJNSdpJMjPPnHM+79drXknPTNtv2uTM+Z7znOcx1loBAAAAALAt83wHAAAAAACEAwUSAAAAAJATCiQAAAAAICcUSAAAAABATiiQAAAAAICclPsOUGp23HFHu3z5ct8xAAAAAMCLtWvXvmStXTbVcxTILSxfvlydnZ2+YwAAAACAF8aYp6Z7jiGsAAAAAICcUCABAAAAADmhQAIAAAAAckKBBAAAAADkhAIJAAAAAMgJBRIAAAAAkBMKJAAAAAAgJxRIAAAAAEBOKJAAAAAAgJxQIAEAAAAAOaFAAgAAAAByQoEEAAAAAOSEAgkAAAAAyAkFEgAAAACQEwokAAAAACAnFEgAAAAAQE4okAAAAACAnFAgAQAAAAA5oUACAAAAAHJCgQQAAAAA5IQCCQAAAADICQUSAAAAAJATCiQAAAAAICcUSAAAAABATiiQAAAAAICcUCABAAAAADmhQAIAAAAAckKBBAAAAADkhAIJAAAAAMgJBRIAAAAAkBMKJAAAAAAgJxRIAAAAAEBOKJAAAAAAgJxQIAEAAAAAOaFAAgAAAAByQoEEAAAAAOSEAgkAAAAAyAkFEgAAAACQEwokAAAAACAnFEgAAAAAQE4okAAAAACAnFAgAQAAAAA5oUACAAAAAHJCgQQAAAAA5IQCCQAAAADICQUSAAAAAJATCiQAAAAAICcUSAAAAABATiiQAAAAAICcUCABAAAAADmhQAIAAAAAckKBBAAAAADkhAIJAAAAAMgJBRIAAAAAkBMKJAAAAAAgJxRIAAAAAEBOKJAAAAAAgJxQIAEAAAAAOaFAAgAAAAByQoEEAAAAAOSEAgkAAAAAyAkFEgAAAACQEwokAAAAACAnFEgAAAAAQE4okAAAAACAnFAgAQAAAAA5oUACAAAAAHJCgQQAAAAA5IQCCQAAAADICQUSAAAAAJATCiQAAAAAICcUSAAAAABATiiQAAAAAICclPsOAAAAnPaujFo7etSbDVSdTKipvkYNdSnfsQAAeA0FEgCAzfgqce1dGTW3dSsYHpUkZbKBmtu6Jalofz/lFQCwLQxhBQBg3ESJy2QDWb1e4tq7MgX/u1s7el4rjxOC4VG1dvQU/O/2+XUDAMKFAgkAwDifJa43G8xoez75/LoBAOFCgQQAYJzPEledTMxoez75/LoBAOFCgQQAYJzPEtdUX6NERdmkbYmKMjXV1xT87/b5dQMAwoUCCQAoOe1dGa1cvUYrzrtZK1evKdq9eD5LXENdSi2NtUolEzKSUsmEWhprizKRjc+v29f/NQBgdpiFFQBQUnzORjrx5/uajbShLuVl5lNfX7fvmWcBADNnrLW+M5SUdDptOzs7fccAgNhauXqNMlPce5dKJvSH8470kAiFwv81AJQmY8xaa216qucYwgoAKClM6BIf/F8DQPhQIAEAJYUJXeKD/2sACB8KJACgpPic0AXFxf81AIQPk+gAAEqK74lsUDz8XwNA+DCJzhaYRAcAAABAnDGJDgAAAABgzhjCCgDYqvauDEMMETl8XwPA7FAgAQDTYqH3mbNWeuUV6bnnpGxW2rhR2rTJPSY+f/VV95pXX5WCQBodlcbGXn8YI5WVSfPmuY/z50tLlkhLl0rJpLRokbR4sfs48fmSJdJOO0nLlrnfh+nxfQ0As0eBBABMq7Wj57WD7AnB8KhaO3pieaA9PCw98YT06KPSs8+6kvjEE9LTT7vPX3zRlcLyclf6Joqcta4kjoy4P2NsbPYZ5s1zf355+eslc+LvGBpyj8WLpR13lHbeWdp1V2nFCimVkqqrpbe+VXrLW6REjFfK4PsaAGaPAgkAmFYcF3ofHZWeesqVxEcflR54QLr/funxx6WXXnLFq6zMFcEgcMVtqj9jcLAw+cbGXi+K03n1Vfd47DH3a2OkBQtcqbVW6u93VzJXrJD22Ud6+9ulPfaQ3vY2t23+/MJkLxVx/L4GgHyhQAIAplWdTCgzxUF1VBZ6Hxx0BfEvf5H++Ef3eOwxV6DKy93zAwOTf8+mTX6yzoW17uvY/Gt5+WX36Ox0X28i4YrvwIC7annQQdI73ym94x3SfvtJVVX+8udb1L+vAaCQKJAAgGk11ddMuldMCu9C78PDUleXdO+90u9/L91zj/TMM9LCha449fdPfm2cbHlF84kn3OPGG6WKCvdvs9NO0gEHSIcdJqXTrmCGdRhslL6vAaDYKJAAgGmFeaH3wUFXFu+4Q7r5Zumvf3VX2kZG3NDTCXErizMRBK//Wz37rHv8+tduOGwQSHvuKR1zjHTUUdKhh4bnKmWYv68BwDdjp7p5I8bS6bTt7Oz0HQMAMEMDA24I6kRhfOABd3UxCCiJhVRW5opjELjJeVatkt79bjf8dbvtfKcDAMyGMWattTY91XNcgQQAhFZvryuL114r3X23uzLW3++uMkpbn2gG+TE6Km3Y4D5/+GHpkUekH/7QFfr99pNOPlk67jg3+ysAIPy4ArkFrkACQOkaG3OTvrS3Sz//ubuHsbxc6uvznQzTWbjQzQK7/fZSY6N7vPOd7t5KAEBp4gokACC0hoel22+XrrpKuuUWVyIHB18fllqo5TKQHxMzvwaBdNll0tVXu/+7I46QPv5x6dhjwzsZDwDEEQUSAEKkvSsTi4k/xsbckNQrrpB+8Qu3bdOmqddcRHhsPtz1lluku+5yw41XrZJOP93dOxnVK5Nx+dkFEH0USAAIifauzKSlBzLZQM1t3ZIUiQNRa6X77nNXqK6+2l1Z7O93ZRLRtHGj+9jWJt16q/se+OAHpVNPlVaulObN85svX6L+swsgXiKyawaA6Gvt6Jm0bp0kBcOjau3o8ZQoP/73f6WvfEV685tdafjWt9wC95s2UR7jZONG939+9dVuaZA3vUn67Gel9et9J5u7qP7sAognCiQAhERvNpjR9lI2NuauOK1aJS1fLl14oVtjcPMZVBFPY2OuSL78svTtb0u1tdIhh7ihzGFdjiVKP7sAQIEEgJCoTk4908h020vRiy9Kq1dLqZSbjfPXv3ZDVScmWgE2Nzzsvjf+/Gc3rHXZMunzn5eefNJ3spmJws8uAEygQAJASDTV1yhRUTZpW6KiTE31NZ4S5cZa6Xe/k044Qdp1V+nLX3bDVjdt8p0MYbJxo/Tqq9Kll0p77eWWArnhBjcxT6kL688uAEyFSXQAICQmJtsIy0yOIyNucpQLLpAyGTc8lVlUMVdDQ+7jH/7gJl1avFj693+XTjmldJcDCdvPLgBsjbG8m0+STqdtZ2en7xgAEFr9/dKVV7qJcfr6uNKIwquqksrKpM98RvrUp6SlS30nAoBwM8astdamp3qOIawAgLx46SV3tXGnnaRzz5Wef57yiOLo63PrS150kbTLLtKZZ0pPPeU7FQBEEwUSADAnzz4rnXGGu7/x4ovdvWp9fb5TIY6CwD1++ENpzz2l979fevBB36kAIFookACAWXnhBemss6Q99pB+9CM3WyazqaIUTMzeetNNUjotnXSS9PjjvlMBQDRQIAEAM5LNuiGqy5e7Kz0DA+Fdnw/RNjrqrkhef720777SP/+zm9AJADB7FEgAQE42bXJLcOyyi/TNb7oD88FB36mAbRsZcd+v//M/0lvfKp19tluTFAAwcxRIAMBWDQ1Jl1ziiuPXv+7ub2SoKsJoYmjrD37grqCfd56bfCdq2rsyWrl6jVacd7NWrl6j9i4uuwLIHwokAORBVA/YbrlF2n13N7vqq6+6JTqAsBscdN/Ll14qvfnNbtmZsTHfqfKjvSuj5rZuZbKBrKRMNlBzW3dk9kkA/KNAAsAcRfGA7ZFHpMMPl0480d0zRnFEFA0MuBMjZ58t1dZKf/6z70Rz19rRo2B4dNK2YHhUrR09nhIBiBoKJADMUZQO2DZskD79aWn//aW77qI4Ih76+txyH0cc4WZsfe4534lmrzcbzGg7AMwUBRIA5igKB2xjY24Y3267ufvDgiA6Q/qAXAWBdMMNbqKdCy8M5yRR1cnEjLYDwExRIAFgjsJ+wPbAA9J++7lhfNmsO4gG4mpoyF15v/BCd//vb3/rO9HMNNXXKFFRNmlboqJMTfU1nhIBiBoKJADMUVgP2IaGpPPPlw46SFq3zg3jA+D09Um9vdL73iedcoq7VzIMGupSammsVSqZkJGUSibU0lirhrqU72gAIsJYa31nKCnpdNp2dnb6jgEgZNq7Mmrt6FFvNlB1MqGm+pqSPmC75x7pH/9ReuEF7nMEtmXhQqmyUrrqKunYY32nAYDCM8astdamp3yOAjkZBRJAlPX1SeeeK11xBUNVgZmqrJTe8x53n/CyZb7TAEDhbK1AMoQVAGJizRo3OQjlEZid/n7pV79yP0c//rHEOXgAcUSBBICICwLptNOk446T/vd/KY/AXAwNueVuzjhDOvpo6eWXfScCgOKiQAJAhHV3S3vvLV17Lfc6AvnU1+dmaH3b26Q77vCdBgCKhwIJABFkrXTppdLBB0tPPcVVR6AQhoakV15xM7Wec440POw7EQAUXskWSGPMk8aYbmPMX40xnePblhpjbjXGPDr+cfvx7cYY801jzHpjzP3GmHds9uecMv76R40xp/j6egCgWF56yU30cf75rjhynxZQWEEgXXaZVFcnPfaY7zQAUFglWyDHHWGt3X+zGYDOk3S7tXYPSbeP/1qSVknaY/xxhqTLJFc4JX1R0sGSDpL0xYnSCQBRdNttbkjdXXexriNQTP390kMPSfvtJ11zje80AFA4pV4gt3SCpKvGP79KUsNm26+2zp8kJY0xO0uql3SrtfYVa+3fJN0q6ehihwaAQhsdlT73OemEE6S//c0NrQNQXGNj7sTNJz8pffCDnMQBEE3lvgNshZX0G2OMlfQ9a+33Jf2dtfY5SbLWPmeMedP4a1OSntns9z47vm267QBipr0ro9aOHvVmA1UnE2qqr1FDXTR2By+/LB1/vPTXvzJRDlAK+vulm2+W9t9f+s1vpBUrfCfKryjvTwFsWykXyJXW2t7xknirMebhrbzWTLHNbmX75N9szBlyQ1/15je/eTZZAZSw9q6Mmtu6FQyPSpIy2UDNbd2SFPqDnvvuk+rrueoIlJqBAenxx92Q1uuuc/clR0GU96cAclOyQ1ittb3jH1+QdL3cPYzPjw9N1fjHF8Zf/qykXTf77btI6t3K9i3/ru9ba9PW2vSyZcvy/aUA8Ky1o+e1g50JwfCoWjt6PCXKj2uvlQ49VHr+ecojUIrGxqSNG93Q8gsvjMaEVlHdnwLIXUkWSGNMlTFm8cTnkt4r6QFJN0qamEn1FEk3jH9+o6R/Gp+N9RBJr44Pde2Q9F5jzPbjk+e8d3wbgBjpzU69hsV020vdyIj0qU9Jp5/OkFUgDILAFcjjjw//fZFR258CmLlSHcL6d5KuN8ZILuOPrbW/NsbcK+lnxpjTJD0t6cTx198i6RhJ6yX1S/pnSbLWvmKM+Yqke8df92Vr7SvF+zIAlILqZEKZKQ5uqpMJD2nm5uWXpWOPle6/n/IIhElfn5sl+e1vd/dFvuUtvhPNTpT2pwBmpySvQFprH7fW7jf+2Mda+7Xx7S9ba4+y1u4x/vGV8e3WWnuWtfYt1tpaa23nZn/WFdbat44/rvT1NQHwp6m+RomKsknbEhVlaqqv8ZRodnp6pH33lf7yF8ojEEYDA9KTT7r1In//e99pZicq+1MAs1eSBRIA8qmhLqWWxlqlkgkZSalkQi2NtaGa8OHuu6WDDuJ+RyDsJu6LrK+Xfv5z32lmLgr7UwBzY2wU7ujOo3Q6bTs7O7f9QgAokrY26aMfdfdRAYiOREL68pelc87xnQQAJjPGrLXWpqd6jiuQAFDC/uu/KI9AVAWB9MUvSmed5a5MAkAYUCABoASNjUlnny2dfz7lEYiy/n7pRz9yM7QODPhOAwDbRoEEgBIzMCA1NEhXXMFkOUAc9PdLa9a4dV1fYa54ACWOAgkAJWTTJunww910/5RHID6CQFq3TnrHO6TeXt9pAGB6FEgAKBHZrLRypXTffQxbBeJoaEjKZKR0WnrqKd9pAGBqFEgAKAEvvywdcoj08MPcBwXE2ciI9MIL0oEHSuvX+04DAG9EgQQAzyYOFp94gjUeAUijo+6k0kEHuZNKAFBKKJAA4NHzz7uDxGefpTwCeN3YmBvWfuih0oMP+k4DAK+jQAKAJ88958pjb680POw7DYBSY+3r90avW+c7DQA4FEgA8OD556WDD6Y8Atg6a6VXX3UlkiuRAEpBue8AALC59q6MWjt61JsNVJ1MqKm+Rg11Kd+x8iqblQ47zF2BHBnxnQZAqbNW2rDB7TfWrpWWL/edqDDisP8HooArkABKRntXRs1t3cpkA1lJmWyg5rZutXdlfEfLm/5+6cgj3RT9lEcAudp8OOvzz/tOk39x2P8DUUGBBFAyWjt6FAyPTtoWDI+qtaPHU6L8GhqSjjlGeughJswBMHNjY9KLL0rvfKcrk1ES9f0/ECUUSAAlozcbzGh7mIyOSiedJN17L+s8Api94WHp6aelI45wIxqiIsr7fyBqKJAASkZ1MjGj7WFhrfSJT0i33hqtAz4AfgwNufUhjzkmOqMZorr/B6KIAgmgZDTV1yhRUTZpW6KiTE31NZ4S5ce550o//SnlEUD+DAxI99wjnXiiG+EQdlHd/wNRRIEEUDIa6lJqaaxVKpmQkZRKJtTSWBvqWfi+/W3pO9+hPALIvyCQbrtNOvts30nmLor7fyCqjLXWd4aSkk6nbWdnp+8YACLg1lulE05wB3kAUCiVldLFF0tnnuk7CYCoMMastdamp3qOK5AAUAAPPSQ1NlIeARRef7/0uc9Jt9/uOwmAOKBAAkCevfSSW+uxr893EgBxEQTS+98v9bDqBYACo0ACQB4NDUlHHy298oqbfRUAimXTJumoo9z+BwAKhQIJAHlirXTKKdKDD0Znan0A4WGt9OKL0qpV7IMAFA4FEgDypKVFuukm7nsE4M/QkNTdLZ12GqMgABQGBRIA8uA3v5G++lXuewTgXxBI11/vlhACgHyjQALAHD37rHTSSVx5BFA6+vqkz39eWrvWdxIAUUOBBIA5GB6WjjuOK48ASk8QSMceK2WzvpMAiBIKJADMwTnnSI88Io2M+E4CAG/0t79J//iP3A8JIH8okAAwSzfeKF1+uVvEGwBK0eCg9Ic/SBdd5DsJgKigQALALDzxhPSRj1AeAZS+vj7pP/7DFUkAmKty3wEAhFt7V0atHT3qzQaqTibUVF+jhrqU71gFNTgove99lEcA4REE0gknSA89JC1b5jtNccTx/QkoBq5AApi19q6Mmtu6lckGspIy2UDNbd1q78r4jlZQn/2s9OST0tiY7yQAkLuNG+NzP2Rc35+AYqBAApi11o4eBcOjk7YFw6Nq7ejxlKjwfvtb6corWbIDQPgMDUn33OPu3Y66OL4/AcVCgQQwa73ZqVvUdNvDbsMGd/ae8gggrPr6pP/7f9193FEWt/cnoJgokABmrTqZmNH2sDvzTFciASDMBgelk06K9jD8uL0/AcVEgQQwa031NUpUlE3alqgoU1N9jadEhfPLX0rt7dLAgO8kADA3o6NuMp1LLvGdpHDi9P4EFBsFEsCsNdSl1NJYq1QyISMplUyopbE2crPcvfSS9LGPMesqgOjo65P+/d9dkYyiuLw/AT4YG4epuGYgnU7bzs5O3zEAlAhrpeOOk2691U1AAQBRYYy0557SffdJFRW+0wAoJcaYtdba9FTPcQUSALbipz91M69SHgFEjbXSU09JX/mK7yQAwoQCCQDTyGbdxDl9fb6TAEBh9PdLF18sPfqo7yQAwoICCQDTOOccJs0BEH2Dg9LHP+6uSALAtlAgAWAKnZ3Sj39MgQQQfWNj7j7In/3MdxIAYUCBBIAtjI66WVcD1psGEBN9fdInPym9+qrvJABKHQUSALbw7W9LzzzjOwUAFNfAgPT5z/tOAaDUUSABYDPPPSedfz4T5wCIn4EB6ZprpL/8xXcSAKWMAgkAmznzTJbsABBfQeCG8I+O+k4CoFRRIAFg3J13SrfeKg0P+04CAP489ZR05ZW+UwAoVRRIAJCbvv7MM92aaAAQZ3197l5IhvIDmAoFEgAk/fzn0tNP+04BAKVhcFC6+GLfKQCUIgokgNgbGpI+9Slp0ybfSQCgNPT3SxddJL3wgu8kAEpNue8AAPxo78qotaNHvdlA1cmEmupr1FCX8h3Li29/m/IIAFsaGZG+8AXp8st9J/GL90tgMq5AAjHU3pVRc1u3MtlAVlImG6i5rVvtXRnf0Youm5W+9CXu9QGALQ0NSf/zP9Ijj/hO4g/vl8AbUSCBGGrt6FEwPHmO9mB4VK0dPZ4S+fOVrzDrKgBMZ3hY+vSnfafwh/dL4I0okEAM9WaDGW2PqmeekS67zC2eDQB4o9FRt8TR3Xf7TuIH75fAG1EggRiqTiZmtD2qvvhFd48PAGB6/f3S5z7nO4UfvF8Cb0SBBGKoqb5GiYqySdsSFWVqqq/xlKj4Mhnp2msZvgoAubj/fulPf/Kdovh4vwTeiAIJxFBDXUotjbVKJRMyklLJhFoaa2M1q9zXviaNjflOAQDhEATSeef5TlF8vF8Cb2Sstb4zlJR0Om07Ozt9xwBQQC+8IC1f7g6IAAC5SSSkP/xBqqvznQRAoRlj1lpr01M9xxVIALHz9a9z9REAZmpw0K0LCSDeKJAAYuWVV6T/9//cgRAAIHdjY9LvfietW+c7CQCfKJAAYuUb3+DqIwDM1tCQdMEFvlMA8IkCCSA2NmxwBZJ1HwFgdkZHpV/9Snr0Ud9JAPhCgQQQG1dcITFvGADMzciIu5ccQDxRIAHEwtiYdNFFbkFsAMDsjYxIP/6xG9UBIH4okABi4Te/kTZu9J0CAKLBGOnKK32nAOADBRJALFx4obRpk+8UABAN/f1uVAeTkgHxQ4EEEHnr10v33us7BQBEy4YNbnQHgHihQAKIvEsucTMHAgDyZ9MmqaXFdwoAxUaBBBBpmzZJV10lDQ/7TgIA0XPPPW6UB4D4oEACiLSrr3aTPQAA8m901K2vCyA+KJAAIu2SS6S+Pt8pACCahofdKI/BQd9JABRLue8AAGanvSuj1o4e9WYDVScTaqqvUUNdynesknL//dJzz/lOAQDRNm+edPPNUmOj7ySlh/dqRBFXIIEQau/KqLmtW5lsICspkw3U3Nat9q6M72gl5fLLpaEh3ykAINo2bpS+8x3fKUoP79WIKgokEEKtHT0KhidPKxoMj6q1o8dTotIzMuLufxwZ8Z0EAKLvD3+QXnzRd4rSwns1oooCCYRQbzaY0fY4uu02FrgGgGIpK5OuvdZ3itLCezWiigIJhFB1MjGj7XF02WVuWBUAoPD6+xnGuiXeqxFVFEgghJrqa5SoKJu0LVFRpqb6Gk+JSsuGDVJHh+8UABAvzzwjPfig7xSlg/dqRBUFEgihhrqUWhprlUomZCSlkgm1NNYys9u4X/xCqqjwnQIA4mVkRLriCt8pSgfv1YgqY631naGkpNNp29nZ6TsGgDlIp6W1a32nAID42WEHN5mOMb6TAJgLY8xaa216que4AgkgUl5+Weru9p0CAOJpcFD6y198pwBQSBRIAJFy883S/Pm+UwBAPA0OutsIAEQXBRJApFxzjbRpk+8UABBPw8Ms5wFEHQUSQGT090t33eU7BQDE2wsvSOvX+04BoFAokAAi47bbpAULfKcAgHizVmpv950CQKFQIAFExrXXujUgAQD+DAy42wkARBMFEkAkjI66CXQAAP719LihrACihwIJIBLuvtt3AgDAhPJy6aabfKcAUAgUSACRcPPNbhIdAIB/fX3Sddf5TgGgECiQACLhppvcMFYAQGn4/e+lsTHfKQDkGwUSQOht3Cg9+qjvFACAzVkrrVvnOwWAfKNAAgi9u+6SEgnfKQAAmxsZkW6/3XcKAPlGgQQQer/+tbRpk+8UAIDNDQxIN9zgOwWAfCv3HQAA5uqWW7jPBgBK0Z//7K5ElnPEOSPtXRm1dvSoNxuoOplQU32NGupSvmMBkiiQQEngjWL2Xn5ZeuYZ3ykAAFMpL5e6uqQDD/SdJDzauzJqbutWMOxmhstkAzW3dUsSxwYoCQxhBTybeKPIZANZvf5G0d6V8R0tFH73O2nhQt8pAABTGRyUbrvNd4pwae3oea08TgiGR9Xa0eMpETAZBRLwjDeKuenocLOwAgBKz9CQW2YJuevNBjPaDhQbBRLwjDeKubnzTjdVPACgNN1/P/vpmahOTj2t+HTbgWKjQAKe8UYxe8PD0mOP+U4BANgaa6UnnvCdIjya6muUqCibtC1RUaam+hpPiYDJKJCAZ7xRzN4DD3D/IwCUurIyqbPTd4rwaKhLqaWxVqlkQkZSKplQS2MtE+igZDALK+DZxBsCs7DOXGenNDq67dcBAPzZtEm6+27ppJN8JwmPhroUxwEoWRRIoATwRjE7d94p9ff7TgEA2Bpr3f4aQDQwhBVAaN19t+8EAIBcPPSQNDbmOwWAfKBAAgilIJCeftp3CgBALsrLpUce8Z0CQD5QIAGE0n33SZWVvlMAAHJ1772+EwDIBwokgFC67z5pZMR3CgBALjZtYiZWICookABCqbubCXQAIEy6unwnAJAPFEgAocSBCACEC/dAAtFAgQQQSuvX+04AAJiJl16ShoZ8pwAwVxRIAKEzOCi9/LLvFACAmaislB57zHcKAHNFgQQQOuvXMwMrAISNMVJPj+8UAOaKAgkgdB5+2DVrqyEAACAASURBVB2IAADCo79feugh3ykAzBUFEkDoPPwwM7ACQNiMjDABGhAFFEgAodPVxRqQABBG69b5TgBgriiQAEKHSRgAIJx6e30nADBX5b4DAMBMPf+87wQAgNnYsEEaHZXKynwniY72roxaO3rUmw1UnUyoqb5GDXUp37EQYRRIoEDYoRfOK6/4TgAAmI0FC6QXXpB23tl3kmho78qoua1bwfCoJCmTDdTc1i1JHHOgYBjCChTAxA49kw1k9foOvb0r4zta6G3cKI2N+U4BAJiNigqGseZTa0fPa+VxQjA8qtYO1ktB4VAggQJgh144mYy0cKHvFACA2cpwLjVverPBjLYD+UCBBAqAHXrh9PZK89hzAUAoDQ9zBTKfqpOJGW0H8oHDMKAA2KEXTm+vm4ABABA+QSA9+6zvFNHRVF+jRMXkGYkSFWVqqq/xlAhxQIEECoAdeuH09kqDg75TAABmi6WY8qehLqWWxlqlkgkZSalkQi2NtUygg4JiFlagACZ23MzCmn9PPumGQAEAwunpp30niJaGuhTHFygqCiRQIOzQC4M1IAEg3LJZ3wkAzAVDWAGEyquv+k4AAJiLjRt9JwAwFxRIAKFCgQSAcOvv950AwFxQIAGEyqZNvhMAAOYiYEUrINQokABChQIJAOE2MCBZ6zsFgNmiQAIIFYY+AUC4GcNyTECYUSABhApDnwAg3CoqmEgHCDMKJIDQsJaz1gAQduXlFEggzCiQAEJjYMANfQIAhNe8eRRIIMwokABCY2TEHXgAAMLLGLc/BxBOkT8UM8YcbYzpMcasN8ac5zsPgNkbHeUKJABEwdiY7wQAZivSBdIYUybpO5JWSdpb0snGmL39pgIwW2NjFEgAiAIKJBBekS6Qkg6StN5a+7i1dkjSTySd4DkTgFmiQAJANFAggfCKeoFMSXpms18/O75tEmPMGcaYTmNM54svvli0cABmhoWnASAaKJBAeJX7DlBgU12reMMhqLX2+5K+L0npdJpDVKBEzZtHiQSAKCgr850gftq7Mmrt6FFvNlB1MqGm+ho11L3hugqwTVEvkM9K2nWzX+8iqddTFkQEO2B/ysookAAQBcyoXVztXRk1t3UrGB6VJGWygZrbuiWJYxjMWNR/fO+VtIcxZoUxZr6kD0m60XMmhNjEDjiTDWT1+g64vSvjO1oscMABANHA/ry4Wjt6XiuPE4LhUbV29HhKhDCL9I+vtXZE0r9K6pD0kKSfWWvX+U2FMGMH7Ne8edw3AwBhZy0Fsth6s8GMtgNbE/UhrLLW3iLpFt85EA3sgP1KJFh8GgDCzlqpstJ3inipTiaUmeJYpTqZ8JAGYcf5H2AGptvRsgMujrIyqaLCdwoAwFyMjEjbbec7Rbw01dcoUTF55qJERZma6ms8JUKYUSCBGWAH7F+Crg4AoTYyIi1e7DtFvDTUpdTSWKtUMiEjKZVMqKWxlgl0MCuRH8IK5NPEjpZZWP1JJKQNG3ynAADM1siIVFXlO0X8NNSlOF5BXlAggRliB+wXBx0AEG7z5zOJDhBm/PgCCBWGPQFAuC1c6DsBgLmgQAIIFQokAIQbM7AC4UaBBBAqS5b4TgAAmAtuRQDCjQIJIFR22MF3AgDAXLCEBxBuFEgAobLbbpIxvlMAAGarutp3AgBzQYEEECq77MIEDAAQZrvv7jsBgLmgQAIIlepqqaLCdwoAwGzMn+9GkgAILwokgFBh6BMAhNeCBezHgbCjQAIIlepqaXjYdwoAwGzMm0eBBMKOAgkgVJYtk4aGfKcAAMzG6CgFEgg7CiSAUCkrYwp4AAirgQFp5519pwAwFxRIAKGzbJnvBACA2Sgvl6qqfKcAMBflvgMAwEztuqv0yCO+UwAAZooTgKWpvSuj1o4e9WYDVScTaqqvUUNdyncslCiuQAIInf32k4zxnQIAMFNve5vvBNhSe1dGzW3dymQDWUmZbKDmtm61d2V8R0OJokACU2jvymjl6jVacd7NWrl6DTvRErPPPlJlpe8UAICZMEaqq/OdAltq7ehRMDw6aVswPKrWjh5PiVDqGMIKbGHiTNzEznTiTJwkhnOUiJoadx8NACA8KivdCUCUlt5sMKPtAFcggS1wJq701dRIAe9rABAq5eVu/43SUp1MzGg7QIEEtsCZuNK3445SRYXvFACAmQgCCmQpaqqvUaKibNK2REWZmur5z8LUKJDAFjgTFw7Ll/tOAACYifnzpaVLfafAlhrqUmpprFUqmZCRlEom1NJYy207mBZ3EQFbaKqvmXQPpMSZuFK0777SunW+UwAAcrVihe8EmE5DXYrCiJxxBRLYAmfiwqGujmGsABAm++7rOwGAfOAKJDAFzsSVvn32kRIJaXjYdxIAwLYsWCCl075TAMgHrkACCKUDDpAGB32nAADkYuFC6cADfacAkA8USAChtPPObk0xAEDp6+93tx4ACD8KJIDQ4mAEAMJhp52kRYt8pwCQDxRIAKF1+OFuYWoAQGk7+GDfCQDkCwUSQGgdfDDDWAGg1CUS7oQfgGigQAIIrXRaGhjwnQIAsDUVFczACkQJBRJAaC1dKiWTvlMAALamv1/abz/fKQDkCwUSQKhxVhsASttuu7llPABEAwUSQKgdfbS7vwYAUHqMkY480ncKAPlEgQQQakceKc1jTwYAJWnxYumYY3ynAJBPHHYBCLW996ZAAkCpGhiQ/uEffKcAkE8cdgEINWM4OAGAUrXbbtL22/tOASCfWIIbQOgdf7x0xx1SX5/vJACACfPmSe97n+8UmIv2roxaO3rUmw1UnUyoqb5GDXUp37HgGQUSscaOMRqOPFIaG/OdAgCwuUWLpPp63ykwW+1dGTW3dSsYHpUkZbKBmtu6JYljpZhjCCtia2LHmMkGsnp9x9jelfEdDTO0++5SVZXvFACAzQWBdNhhvlNgtlo7el4rjxOC4VG1dvR4SoRSQYFEbLFjjA5jpKOO8p0CALC5Pffk5F6Y9WaDGW1HfFAgEVvsGKPluOPcdPEAAP8qKtz96Qiv6uTUiyxPtx3xQYFEbLFjjJZVq6ShId8pAACStGCB1NjoOwXmoqm+RomKsknbEhVlaqqv8ZQIpYICidhixxgtS5dKtbW+UwAAJFcg6+p8p8BcNNSl1NJYq1QyISMplUyopbGWCXTALKyIr4kdILOwRsdHPyqtW+cmbgAA+DFvnrv6aIzvJJirhroUx0V4A2Ot9Z2hpKTTadvZ2ek7BoBZeOopN2nDwIDvJAAQX9ttJ7W1MbkZEGbGmLXW2vRUzzGEFUBk7LablOJEKQB4NTIivetdvlMAKBQKJIBI+chHpPnzfacAgPiqr3ezsAKIJgokgEj5wAc4cAEAXxYvdifyAEQXBRJApNTWsnA1APgyNOSuQAKILgokgEgxxs3GylVIACi+d71LWrTIdwoAhUSBBBA5p54qlbNIEQAU1eLF0lln+U4BoNAokAAiZ599pF139Z0CAOLFWmnVKt8pABQaBRJAJJ11llRZ6TsFAMRDWZn04Q8zCzYQBxRIAJF08snS6KjvFAAQDwsXSp/4hO8UAIqBAgkgkpYtkw491HcKAIiH7beXDjjAdwoAxUCBBBBZZ53lJnUAABTOwoXSJz/pZsEGEH0USACRdeyxDGMFgGL4p3/ynQBAsVAgAUTWggVuUgeW9ACAwqmrY+brOGnvymjl6jVacd7NWrl6jdq7Mr4jocgokAAi7XOfkyoqfKcAgGhavFhqbvadAsXS3pVRc1u3MtlAVlImG6i5rZsSGTMUSMQKZ83iZ889pf32850CAKIpkZCOOcZ3ChRLa0ePguHJ94YEw6Nq7ejxlAg+UCARG5w1i68vfEFatMh3CgCIlspK6Zxz3BqQiIfebDCj7YgmCiRig7Nm8XXMMe5ABwCQP2Nj0umn+06BYqpOJma0HdFEgURscNYsvsrK3L2QlEgAyI+yMunEE936j4iPpvoaJSomX3JOVJSpqb7GUyL4QIFEbHDWLN5OP92dLQcAzN38+VJTk+8UKLaGupRaGmuVSiZkJKWSCbU01qqhLuU7GoqIye0RG031NWpu6540jJWzZvGxdKn0wQ9K117L2pAAMFd77SXV1vpOAR8a6lIUxpjjCiRig7NmOPdcd9YcADB7ixZJ//ZvvlMA8MVYa31nKCnpdNp2dnb6jgGgQN79bmnNGoldHwDMzpvfLD3xhDSPyxBAZBlj1lpr01M9x48+gFi58EK3bhkAYOYWLZJaWiiPQJzx4w8gVg46SNp/f98pACCcFi+WTjrJdwoAPlEgAcTO6tVSVZXvFAAQLosWSV/5ilTOFIxArFEgAcTOYYdJNUy+CwAzsnCh9LGP+U4BwDcKJIBY+vrXuQoJALmqqpK+9CVmsgZAgQQQU0cdJS1f7jsFAIRDRYV02mm+UwAoBRRIALFkjHTRRVyFBIBtqayULrjADWEFAAokgNhatUraay9XJgEAU1u0SDrrLN8pAJQKCiSA2DJG+u53OasOANOpqpIuvVRasMB3EgClggIJINYOPFB6z3uYlh4AprLrrqz7CGAyCiSA2Puv/3ITRAAAXldZKV12mTSPo0UAm2GXACD2VqyQTj2VoawAMGHePOngg6XDD/edBECpoUACgKT/+A+prMx3CgAoDQsWSN/6lu8UKHXtXRmtXL1GK867WStXr1F7V8Z3JBQBBRKRx84NudhhB+n8892QLQCIs4oKqbFR2mcf30lQytq7Mmpu61YmG8hKymQDNbd1c5wVAxRIRBo7N8zEZz4jbbed7xQA4FdFhdTa6jsFSl1rR4+C4dFJ24LhUbV29HhKhGKhQCLS2LlhJhYulK68kquQAOKrqkr62teknXf2nQSlrjcbzGg7ooMCiUhj54aZOvpo6YgjmJUVQDztsot09tm+UyAMqpOJGW1HdFAgEWns3DAb3/ueNH++7xQAUFyJhHTNNUwohtw01dcoUTH5myVRUaam+hpPiVAsFEhEGjs3zEYq5WZlrarynQQAimPBAulDH5IOPNB3EoRFQ11KLY21SiUTMpJSyYRaGmvVUJfyHQ0FZqy1vjOUlHQ6bTs7O33HQB61d2XU2tGj3myg6mRCTfU17NywTSMj0t57S48+6jsJABTekiXSk09KyaTvJABKgTFmrbU2PdVz5cUOAxRbQ12KwogZKy+Xrr5aOvJIKeCWWQARVlXl1nykPALIBUNYAWAahxwinXSSm50VAKKorMyt9/jRj/pOAiAsKJAAsBWXXiotWuQ7BQAUxoIF0rXXSsb4TgIgLCiQALAVS5ZIP/0pa0MCiJ6qKuk//1PafXffSQCECQUSALbhyCOlj33MTXEPAFFQXi4dcID0L//iOwmAsKFAAkAOLrlEWrrUdwoAyI9EgqGrAGZnqwXSGPP3xpjvGGPuN8a8aIx52hhzizHmLGPMkmKFBADfKiul667jKiSA8KuslL73Pam62ncSAGE0bYE0xvxK0umSOiQdLWlnSXtL+jdJCyXdYIw5vhghAaAUHHywdPbZ3A8JILzmz3fD8k8+2XcSAGFlrLVTP2HMjtbal7b6m3N4Tdik02nb2dnpOwaAEjU0JNXWSo8+Kk2z+wSAkpVMSuvXSzvs4DsJgFJmjFlrrU1P9dy0VyC3LIbGmO2MMUsnHlO9BgCibv586frrGcoKIHwSCeknP6E8ApibbU6iY4z5F2PM85Lul7R2/MElOgCxtffe0ne/y1BWAOFRWemG4NfX+04CIOxymYX1HEn7WGuXW2tXjD9YMQhArJ1yivSBD3AlEkDpq6iQ9t1XuvBC30kAREF5Dq95TFJ/oYMAQNh873vSH/8oPf64NDbmOw0ATG3RIumGG6SyMt9JEEXtXRm1dvSoNxuoOplQU32NGupSvmOhgHIpkM2S7jbG/FnS4MRGa+2nCpYKAEIgkZBuuUWqq5P6+nynAYA3SiTcfds77eQ7CaKovSuj5rZuBcOjkqRMNlBzW7ckUSIjLJchrN+TtEbSn/T6PZBrCxkKyIf2roxWrl6jFefdrJWr16i9K+M7EiJojz2kK6/kfkgApaeqSmpulv7hH3wnQVS1dvS8Vh4nBMOjau3o8ZQIxZDLFcgRa+1nC54EyCPOiKGYTjxRuv126ZprpH4G/AMoAfPnSwceKJ1/vu8kiLLebDCj7YiGXK5A3mGMOcMYs/OWy3gApYozYii2b37TXY2sqPCdBEDcGSMtXSr94hfSvFyO9IBZqk5OPZPcdNsRDbnsVj6s8fsgxTIeCAnOiKHY5s+XfvMbd9BmjO80AOKsqkpas4b1HlF4TfU1SlRMnp0pUVGmpvoaT4lQDLkMYd3dWms332CMWVigPEBeVCcTykxRFjkjhkJ605vcUNaDD2ZSHQB+JBLSdddJe+3lOwniYOK2IGZhjZdcCuQPJZ068QtjTJWkGyUdVahQwFw11ddMugdS4owYimOffaSf/9ytERlwwRtAEVVWSi0t0nvf6zsJ4qShLkVhjJlchrBmjDGXSZIxZntJt0r674KmAuaooS6llsZapZIJGUmpZEItjbXs4FAUq1ZJX/2qG0YGAMVQWSl9+MPSp1hkDUCBmS1Gp079ImO+LmmJpAMkrbbWXlfoYL6k02nb2cktngDmxlrp1FOln/2MmVkBFNaCBW7G1TvukMpzGVsGANtgjFlrrU1P9dy0VyCNMY0TD0n3SDpEUpckO74NADANY6Tvf196+9vdBDsAUAjz5kl/93fSL39JeQRQHFvb1Ry3xa+7JFWMb7eS2goVCgCioKJCuuUWKZ2Wnn5aGhnxnQhA1CxZ4q48LlniOwmAuJi2QFpr/7mYQQAgirbfXvr976UDDpCef14aG/OdCEBULF4s3XmntPvuvpMAiJOtDWH9t/FJc6Z7/khjzLGFiQUA0bHzzq5EJpOsEQkgP6qq3Nqz++7rOwmAuNnaENZuSb80xgxI+oukFyUtlLSHpP0l3SbpwoInBIAI2H136Xe/kw49VNq40XcaAGGWSEjXXy8dcojvJADiaNorkNbaG6y1KyV9UtI6SWWSNsgt4XGQtfYz1toXixMTAMJv333dFQOW9wAwW4mEdNVV0nve4zsJgLja5nxd1tpHJT1ahCwAEHmHHCK1t0vHHy8Fge80AMKkslL6xjekE0/0nQRAnE17BRIAUBjvfrd0zTXuSgIA5KKyUjr/fOmMM3wnARB3FEgA8OADH5Auv5wSCWDbKiulpibpC1/wnQQAchjCCgAojA9/2C0CfuqpDGcFMLXKSum886QLLvCdBACcbRZIY8wySZ+QtHzz11trTy1cLACIhw99yJXIj3+cEglgsolhq1x5BFBKcrkCeYOku+SW7RgtbBwAiJ+TTpIqKqSPfIQSCcBJJKQvfckNXQVKXXtXRq0dPerNBqpOJtRUX6OGupTvWCiQXApkpbX23IInAfKMnRnC5P3vl667TvrgB6X+ft9pAPiUSEgXXyz9n//jOwmwbe1dGTW3dSsYdteZMtlAzW3dksRxV0TlMonOL40xxxQ8CZBHEzuzTDaQ1es7s/aujO9owLRWrZJuvtkNWwMQT4mEdNlllEeER2tHz2vlcUIwPKrWjh5PiVBouRTIT8uVyAFjzMbxx4ZCBwPmgp0Zwurww6U1a6TttpOM8Z0GQDElEtJ//7d0yim+kwC5681Ofe/FdNsRftsskNbaxdbaedbaheOfL7bWbleMcMBssTNDmB18sNTZKe20k7s3EkC0GeNOGt12m9TY6DsNMDPVyanXo5puO8Ivp3UgjTHHG2MuHn8cW+hQwFyxM0PY7bGH9Ne/uo8LF/pOA6BQysulN71Juuce6dBDfacBZq6pvkaJirJJ2xIVZWqqr/GUCIW2zQJpjFktN4z1wfHHp8e3ASWLnRmiYPODSu6LBKJn4ULprW+V7rtPquHtCSHVUJdSS2OtUsmEjKRUMqGWxlom0IkwY63d+guMuV/S/tbasfFfl0nqsta+vQj5ii6dTtvOzk7fMZAHzMKKqBgZkU4/XfrFL6S+Pt9pAORDZaV00EHSTTdJixb5TgMAkxlj1lpr01M9l8syHpKUlPTK+OdL8pIKKLCGuhSFEZFQXi5deaX0lrdIq1ezzAcQdpWV0gc+IF1xhfv5BoAwyeUeyBZJXcaYHxljrpK0VtKFhY0FANicMdIFF0iXX85wViDMEgmpuVm66irKI4Bw2uauy1p7rTHmt5IOlGQknWut/d9CBwMAvNHJJ0t77SXV10vZrDQ05DsRgFzMmydVVbmh6O99r+80ADB7016BNMbsOf7xHZJ2lvSspGckVY9vKwhjzJeMMRljzF/HH8ds9lyzMWa9MabHGFO/2fajx7etN8act9n2FcaYPxtjHjXG/NQYM79QuQGgWPbfX3rwQemAA7gaCYTBwoXSihVuZmXKI4Cw29oQ1s+Of/zPKR4XFzjXN6y1+48/bpEkY8zekj4kaR9JR0v6rjGmbHxSn+9IWiVpb0knj79Wkr4+/mftIelvkk4rcG4AKIoddpDuvNNNrkOJBEpXVZUbMXDffdLuu/tOAwBzN+0QVmvtGeOfrrLWDmz+nDHGx6pkJ0j6ibV2UNITxpj1kg4af269tfbx8Ww/kXSCMeYhSUdK+vD4a66S9CVJlxU1NQAUSHm5dOml0t//vXTaaUyuA5SaRMLdu/z5z7v7mAEgCnKZROfuHLfl078aY+43xlxhjNl+fFtKbgjthGfHt023fQdJWWvtyBbb38AYc4YxptMY0/niiy/m8+sAgIL70IekP/5R2mknacEC32kAlJdL220n3XijdO65lEcA0bK1eyB3MsYcIClhjKkzxrxj/HG4pDkNmDLG3GaMeWCKxwlyVwjfIml/Sc/JDZmV3AQ+W7Kz2P7GjdZ+31qbttamly1bNuOvBwB8e/vbpYceko4+miGtgE9VVe4+5XXrpHe/23caAMi/rc3CWi/p45J2kXTJZts3SvrCXP5Sa21Ou1RjzA8k/XL8l89K2nWzp3eR1Dv++VTbX5KUNMaUj1+F3Pz1ABA5yaR0/fVuzcizz5YGBqSxMd+pgPhIJNxw1fPPl8rKfKcBgMLY2j2QV0m6yhjzAWvtdcUKZIzZ2Vr73Pgv3y/pgfHPb5T0Y2PMJZKqJe0h6R65K417GGNWSMrITbTzYWutNcbcIemDkn4i6RRJNxTr6wAAH4yRTj1VOuww6fjjpaef5t5IoNAWLHAncNrbpUMO8Z0GAApr2gJpjPmotfa/JS03xnx2y+ettZdM8dvy4SJjzP5yw02flPQv43/fOmPMzyQ9KGlE0lnW2tHxrP8qqUNSmaQrrLXrxv+scyX9xBjzVUldkn5YoMwAUFL22MPN+tjUJP3gB1IQ+E4ERFNlpZtl9corpSVLfKcBgMLb2hDWqvGPi4oRZIK19mNbee5rkr42xfZbJN0yxfbH9fpMrQAQK/Pnu1lajztOOukkqa9PGhrynQqIhrIyt77jd74jnXKK7zSAP+1dGbV29Kg3G6g6mVBTfY0a6qactxIRYaydcl6Z2Eqn07azs9N3DADIq5decmtG3norQ1qBuaqqkvbZR7r2WtZ2RLy1d2XU3NatYHj0tW2JijK1NNZSIkPOGLPWWpue6rltLuNhjLnIGLOdMabCGHO7MeYlY8xH8x8TyL/2roxWrl6jFefdrJWr16i9K+M7EuDFjju6+7N+8hNp6VJ35QTAzJSXu/J4ySXSn/5EeQRaO3omlUdJCoZH1drR4ykRiiGXdSDfa63dIOlYuZlQ3yapqaCpgDyYOCuWyQaykjLZQM1t3ZRIxNpxx0mPP+6GtLLcB5C7qirpXe+SenqkM85gbUdAknqzU99gP912REMuBbJi/OMxkq611r5SwDxA3nBWDJjakiXSVVdJv/yltNNObukBAFObP1/abjs3GdVtt0kpRuUBr6lOTv0GMt12REMuBfImY8zDktKSbjfGLJM0UNhYwNxxVgzYuiOOkB57TDrtNFciuaICTFZZKR1zjPs5OflkfkaALTXV1yhRMXnR00RFmZrqazwlQjFss0Baa8+T9PeS0tbaYUl9kk4odDBgrjgrBmxbZaX0rW9Jd90l1da6YXpA3FVVSbvuKl1/vXvsuKPvREBpaqhLqaWxVqlkQkZSKplgAp0Y2NoyHpIkY0zF/2/v7qMrS+s60X9/XZ2mT3ejkXcqiCA0WQusiwUlL4MvYysGX9ZQ1MBSxzWilyvCwDg63mjXONdRQau1ZnSu3vGFARY4o4OIReBKQ4mgMII4NAYtEOp2IyidaoRWgrwEqK5+7h/npDupTlXtpJI65+R8PmudleQ5+5w86d519v4++9m/J8m/TPL11R96e1uSX9/hfsFFm5+b3bAymFExuKfHPz5573uT3/7t5F//6+Tzn7d2JJPniiuSqankp386+aEf6n8PnN/B/TMC44TpMoX115I8PsmvDh6PG7TBSDMqBptTlXzP9yR/+7fJC17Qr9S6Z8+FXwfjrqo/jftZz+oXmfrRHxUeAc7lgutAVtVftNYee6G23cI6kAB9H/pQ8rznJe98p7Uj2b2uvjp5xCOSl7+8fyUegItcBzLJmap6xJo3+8okZ86zPQC7wCMekbz5zf31Ix/2MPdHsrtcfXVy3/smv/Eb/enbwiNAN10C5HySP6qqP66qtyV5a5If3dluATAqnvrU5JZbkpe8JHnwgwVJxttVV/WXsnnxi5Olpf60bdVVAbq74BTWJKmqeyWZTVJJPtha+8JOd2xYTGEFOLfTp5NXvCI5fLhfaOeznx12j6CbXq9/T+/11yc//MMGQgDO56KmsFbVlUlekOSnkvxkkucP2gCYMFNTyQ/8QHLrrcmLXtS/knPVVcPuFZzbve7V30f/zb9JPvrR5Cd+QngEuBhdprD+ZpLHJPmVJP9Pkkcn+W872SkARtuVVyY/8iP9KYD/7t8l11wjVNcx7AAAIABJREFUSDJarryy//iBH0g+8pHkyJFkenrYvQIYf10C5Gxr7TmttT8aPJ6b5FE73TEARt/VV/ev6Nx2W/KzP9tfcP2aa4bdKybZ1Vf398Ef/dH+kjS/8ivJ/e8/7F4B7B5dAuRiVT1p9YeqemKSd+xclwAYN9dc07+v7NSpflXLRzxCkOTSuuaa5AEPSG64Ifm7v+sXyREcAbZflwD5xCTvrKqPVNVHkvxpkm+oqhNV9Zc72jsAxsrUVPIv/kVy883J7/1e8sQn9qe2XtblaANbcPXVybXXJi99aX9K9QtfaDo1wE66vMM2T9vxXgCwq1Ql3/It/cef/3n/atAb39hvX1kZdu8Yd1dc0R+UeMITkp/8yeS66yzFAXCpXDBAttb+5lJ0BIDd6XGPS44dS26/PXn5y5P//J+TT386+cxnht0zxs29791fiuMHfzB5/vOTr/iKYfcIJtfC4lKOHj+ZU8sr2Tvdy/zcbA7unxl2t7gETCoC4JK43/2SH/ux/jTD170u+fZv7y+xcKWFoTiPqan+PvLkJyevfGXy8Y/373MUHmF4FhaXcvjYiSwtr6QlWVpeyeFjJ7KwuDTsrnEJdJnCCmPPKBmMjqr+lMPrrusXO3npS5Nf/uX+1NbPfja5885h95BR8CVf0v/6nOckL3hBvzATMBqOHj+ZldNn1rWtnD6To8dPOr+aABe8AllVL6yqL7sUnYGdYJQMRtcDH3j3MiA33ph83/f1pyne+97D7hnDcO97J71ecuhQ8qpX9ac9/+IvCo8wak4tb3wz+7na2V26TGF9UJJ3V9Wrq+ppVW5TZ7ycb5QMGA2XXZZ87dcmL3tZ8vd/n7z61ckzn9kPE8Lk7nbNNf2pzE99av///+239yv4fuu39qevAqNn73RvU+3sLhcMkK21f5/k2iQvS/J9SW6uqp+rKuOBjAWjZDBepqaSpz0t+d3f7YfJl7+8X831Xve6e1oj4+3e9+5XUn3iE/vTlz/2seQP/iB51rMswQHjYH5uNr2pPevaelN7Mj83O6QecSl1ugeytdaq6mNJPpbkjiRfluQ1VfXm1tqP7WQH4WLtne5laYOwaJQMRl+v178S+cxn9qu2vuUt/WD5+7+fnDmTfOELyenTw+4lF7JnT3+9xjvuSL75m5Pv/M5kbi65732H3TNgK1bvc1RfYjJVa+38G1T9UJJnJ7k9yUuTLLTWTlfVZUlubq3tqiuRBw4caDfddNOwu8E2Wr0Hcu001t7Unhw5tM8HHYypO+9MFheT178++Z3fST7ykeTyy/tFeBgNvcEY3f3ul/zzf5484xnJP/kn/f9PAIy2qnpPa+3ARs91+Ri/X5JDZ68H2Vq7s6q+Yzs6CDvJKBnsPpddljz+8f3HT/90fwrkjTcmr31t8id/knzxi/1qrwLlpdPr9a80ttafmvr0pyff8R3JV37lsHsGwHa64BXISeMKJMB4ay350IeSt70teeMb+18/+9l+uPnMZ4bdu93jqqv6If2yy5KnPKW/ruc3fEPymMf02wAYXxd7BRIAxkZV8shH9h/PeU6/7W/+Jnn725M3vSl55zuTW2/tB6A77kg+97nh9nccXHllv+jNykp/SuoTnpB827clX//1yexs/785AJPBFcizuAIJsPt9/vPJiRP9+yjf8Y7kXe9KPvzhflBqrX/FclIPj1dd1b9P8XOfSx7ykOTAgeTrvi553OOSxz7WsioAk+B8VyAFyLMIkACT6fTp5IMfTP78z5P3vS/5i79Ibr45WVrqB6orruhvsxuuWF55ZX9ZlDvu6FeyfeAD+1dsH/vY/hTU/fuTffv62wEweUxhBYALmJrqh6Z9+9a333lncupUP0zefHPyV3/VD5e33tpf9P4f/7EfLq+4oj+V8447+lc4z5zZ+PfspMsu64e+qan+FdTTp/sB8ZprkvvcJ5mZSb7qq/qPRz0qufba5KEP7d8fCgBdCJAAcB6XXdafyvmQhyTf+I33fP7OO/tB8rbb7n587GP9pUVuvz351Kf6IfMzn+lPjf3sZ/sB8/Of74e8yy7rB8/Vx6rVCUKt3R1Gr7yy/7jqqv66itdc059S+qVf2l9T8Su+Itm7N3nwg5MHPaj/9QEPsHQGANvHIQUALsJll/VD2gMe0J8Cuhlf/GL/ceed/aB45539sHjZZesfU1N3X+EEgGESIAFgSFanvgLAuLBSEwAAAJ0IkAAAAHRiCisAAHBBC4tLOXr8ZE4tr2TvdC/zc7M5uH9m2N3iEhMgAQCA81pYXMrhYyeycrpfFnppeSWHj51IEiFywgiQTByjZwAAm3P0+Mm7wuOqldNncvT4SedRE0aAZKIYPQMA2LxTyyubamf3UkSHiXK+0TMAADa2d7q3qXZ2LwGSiWL0DABg8+bnZtOb2rOurTe1J/Nzs0PqEcMiQDJRjJ4BAGzewf0zOXJoX2ame6kkM9O9HDm0zy1AE8g9kEyU+bnZdfdAJkbPAAC6OLh/RmBEgGSyrH7oqcIKAACbJ0AycYyeAQDA1rgHEgAAgE4ESAAAADoRIAEAAOhEgAQAAKATARIAAIBOBEgAAAA6ESABAADoRIAEAACgk8uH3QEAAGD0LCwu5ejxkzm1vJK9073Mz83m4P6ZYXeLIRMgAQCAdRYWl3L42ImsnD6TJFlaXsnhYyeSRIiccAIkE8/oGgDAekePn7wrPK5aOX0mR4+fdJ404QRIJprRNQCAezq1vLKpdiaHIjpMtPONrgEATKq9071NtTM5BEgmmtE1AIB7mp+bTW9qz7q23tSezM/NDqlHjAoBkolmdA0A4J4O7p/JkUP7MjPdSyWZme7lyKF9bvHBPZBMtvm52XX3QCZG1wAAkn6IFBg5mwDJRFv9UFSFFQAALkyAZOIZXQMAgG7cAwkAAEAnAiQAAACdCJAAAAB0IkACAADQiQAJAABAJwIkAAAAnQiQAAAAdCJAAgAA0Mnlw+4AAAAwXAuLSzl6/GROLa9k73Qv83OzObh/ZtjdYgQJkLABH6IAwKRYWFzK4WMnsnL6TJJkaXklh4+dSBLnP9yDKaxwltUP0aXllbTc/SG6sLg07K4BAGy7o8dP3hUeV62cPpOjx08OqUeMMgESzuJDFACYJKeWVzbVzmQTIOEsPkQBgEmyd7q3qXYmmwAJZ/EhCgBMkvm52fSm9qxr603tyfzc7JB6xCgTIOEsPkQBgElycP9Mjhzal5npXirJzHQvRw7tU0CHDanCCmdZ/bBUhRUAmBQH988416ETARI24EMUAADuyRRWAAAAOhEgAQAA6ESABAAAoBMBEgAAgE4ESAAAADoRIAEAAOhEgAQAAKAT60ACAMAEWVhcytHjJ3NqeSV7p3uZn5u1/jWdCZAAADAhFhaXcvjYiaycPpMkWVpeyeFjJ5JEiKQTARI6MloHAIy7o8dP3hUeV62cPpOjx086r6ETARI6MFoHAOwGp5ZXNtUOZ1NEBzo432gdAMC42Dvd21Q7nE2AhA6M1gEAu8H83Gx6U3vWtfWm9mR+bnZIPWLcCJDQgdE6AGA3OLh/JkcO7cvMdC+VZGa6lyOH9rklh87cAwkdzM/NrrsHMjFaBwCMp4P7ZwRGtkyAhA5WP2RVYQUAYJIJkNCR0ToAACadeyABAADoRIAEAACgEwESAACATgRIAAAAOhEgAQAA6ESABAAAoBMBEgAAgE6sAwkAALvQwuJSjh4/mVPLK9k73cv83Kw1rbloAiRsAx/QAMAoWVhcyuFjJ7Jy+kySZGl5JYePnUgS5yhcFFNY4SKtfkAvLa+k5e4P6IXFpWF3DQCYUEePn7wrPK5aOX0mR4+fHFKP2C0ESLhIPqABgFFzanllU+3QlQAJF8kHNAAwavZO9zbVDl0JkHCRfEADAKNmfm42vak969p6U3syPzc7pB6xWwiQcJF8QAMAo+bg/pkcObQvM9O9VJKZ6V6OHNqngA4XTRVWuEirH8SqsAIAo+Tg/hnnI2w7ARK2gQ9oAAAmgSmsAAAAdCJAAgAA0IkACQAAQCcCJAAAAJ0IkAAAAHQiQAIAANCJZTwAAGBMLSwuWYuaS0qABACAMbSwuJTDx05k5fSZJMnS8koOHzuRJEIkO0aAhB1mZBAA2AlHj5+8KzyuWjl9JkePn3SuwY4RIGEHGRkEAHbKqeWVTbXDdlBEB3bQ+UYGAQAuxt7p3qbaYTsIkLCDjAwCADtlfm42vak969p6U3syPzc7pB4xCQRI2EFGBgGAnXJw/0yOHNqXmeleKsnMdC9HDu1zmww7yj2QsIPm52bX3QOZGBkEALbPwf0zAiOXlAAJO2j1A10VVgAAdgMBEnaYkUEAAHYL90ACAADQyVACZFU9q6reX1V3VtWBs547XFW3VNXJqppb0/60QdstVXX9mvaHV9WfVdXNVfU7VXXFoP1eg59vGTz/sEv19wEAAOxGw7oC+b4kh5K8fW1jVT06yXcleUySpyX51araU1V7kvyXJN+a5NFJvnuwbZL8fJJfaq1dm+STSZ4zaH9Okk+21h6Z5JcG2wEAALBFQwmQrbUPtNY2Wkn96Ule1Vr7Qmvtw0luSfKEweOW1tpft9a+mORVSZ5eVZXkuiSvGbz+lUkOrnmvVw6+f02SbxpsDwAAwBaM2j2QM0k+uubnWwdt52q/b5Ll1todZ7Wve6/B858abH8PVfXcqrqpqm76xCc+sU1/CgAAwO6yY1VYq+oPkzxog6d+orX2unO9bIO2lo2DbjvP9ud7r3s2tvaSJC9JkgMHDmy4DQAAXGoLi0uWA2Ok7FiAbK198xZedmuSL1/z80OSnBp8v1H77Ummq+rywVXGtduvvtetVXV5ki9N8g9b6BPsKAcGAGAjC4tLOXzsRFZOn0mSLC2v5PCxE0niXIGhGbUprK9P8l2DCqoPT3Jtkv+V5N1Jrh1UXL0i/UI7r2+ttSR/lOSZg9c/O8nr1rzXswffPzPJWwfbw8hYPTAsLa+k5e4Dw8Li0rC7BgAM2dHjJ+8Kj6tWTp/J0eMblRKBS2NYy3g8o6puTfLkJG+oquNJ0lp7f5JXJ/mrJG9K8oLW2pnB1cUXJjme5ANJXj3YNkl+PMm/rapb0r/H8WWD9pclue+g/d8muWvpDxgVDgwAwLmcWl7ZVDtcCjs2hfV8WmuvTfLaczz3s0l+doP2G5PcuEH7X6dfpfXs9s8nedZFdxZ2kAMDAHAue6d7WdrgnGDvdG8IvYG+UZvCChPlXAcABwYAYH5uNr2pPevaelN7Mj83O6QegQAJQ+XAAACcy8H9MzlyaF9mpnupJDPTvRw5tE8BHYZqKFNYgb7VA4AqrADARg7un3FewEgRIGHIHBgAABgXprACAADQiQAJAABAJwIkAAAAnQiQAAAAdCJAAgAA0IkqrAAAMCQLi0uW82KsCJAAADAEC4tLOXzsRFZOn0mSLC2v5PCxE0kiRDKyBEgYcUYmAWB3Onr85F3hcdXK6TM5evykYz0jS4CEEWZkEgB2r1PLK5tqh1GgiA6MsPONTAIA423vdG9T7TAKBEgYYUYmAWD3mp+bTW9qz7q23tSezM/NDqlHcGECJIwwI5MAsHsd3D+TI4f2ZWa6l0oyM93LkUP73KbCSHMPJIyw+bnZdfdAJkYmAWA3Obh/RmBkrAiQMMJWDyiqsAIAMAoESBhxRiYBABgV7oEEAACgEwESAACATgRIAAAAOhEgAQAA6EQRHQAA2AELi0sqqbPrCJCwCzhAAcBoWVhcWreW89LySg4fO5EkjtGMNVNYYcytHqCWllfScvcBamFxadhdA4CJdfT4ybvC46qV02dy9PjJIfUItocACWPOAQoARs+p5ZVNtcO4ECBhzDlAAcDo2Tvd21Q7jAsBEsacAxQAjJ75udn0pvasa+tN7cn83OyQegTbQ4CEMecABQCj5+D+mRw5tC8z071UkpnpXo4c2qeADmNPFVYYc6sHIlVYAWC0HNw/43jMriNAwi7gAAUAwKVgCisAAACdCJAAAAB0IkACAADQiQAJAABAJ4roAADAJi0sLqmAzkQSIGFCONABwPZYWFzK4WMnsnL6TJJkaXklh4+dSBLHVnY9U1hhAqwe6JaWV9Jy94FuYXFp2F0DgLFz9PjJu8LjqpXTZ3L0+Mkh9QguHQESJoADHQBsn1PLK5tqh91EgIQJ4EAHANtn73RvU+2wmwiQMAEc6ABg+8zPzaY3tWddW29qT+bnZofUI7h0BEiYAA50ALB9Du6fyZFD+zIz3UslmZnu5cihfQroMBFUYYUJsHpAU4UVALbHwf0zjqNMJAESJoQDHQAAF8sUVgAAADoRIAEAAOhEgAQAAKAT90ACAMDAwuKSonNwHgIksI4DJwCTamFxKYePncjK6TNJkqXllRw+diJJHAthwBRW4C6rB86l5ZW03H3gXFhcGnbXAGDHHT1+8q7wuGrl9JkcPX5ySD2C0SNAAndx4ARgkp1aXtlUO0wiARK4iwMnAJNs73RvU+0wiQRI4C4OnABMsvm52fSm9qxr603tyfzc7JB6BKNHgATu4sAJwCQ7uH8mRw7ty8x0L5VkZrqXI4f2KaADa6jCCtxl9QCpCisAk+rg/hnHPTgPARJYx4ETAIBzMYUVAACATgRIAAAAOjGFFQCAXW9hcck9/rANBEgAAHa1hcWlHD52IiunzyRJlpZXcvjYiSQRImGTBEhgS4zkAjAujh4/eVd4XLVy+kyOHj/p2AWbJEACm2YkF4Bxcmp5ZVPtwLkpogNs2vlGcgFg1Oyd7m2qHTg3ARLYNCO5AIyT+bnZ9Kb2rGvrTe3J/NzskHoE40uABDbNSC4A4+Tg/pkcObQvM9O9VJKZ6V6OHNrntgvYAvdAAps2Pze77h7IxEguAKPt4P4ZgRG2gQAJbNrqAVgVVgCAySJAAltiJBcAYPK4BxIAAIBOXIEEAGAsLSwuuZ0CLjEBEgCAsbOwuLSuoNvS8koOHzuRJEIk7CABEthxRogB2G5Hj59cVw08SVZOn8nR4ycdY2AHCZDAjjJCDMBOOLW8sql2YHsoogPsqPONEAPAVu2d7m2qHdgeAiSwo4wQA7AT5udm05vas66tN7Un83OzQ+oRTAYBEthRRogB2AkH98/kyKF9mZnupZLMTPdy5NA+t0fADnMPJLCj5udm190DmRghBmB7HNw/IzDCJSZAAjtq9cCuCisAwPgTIIEdZ4QYAGB3ECABABgZ1g6G0SZAAgAwEqwdDKNPgARGipFngMl1vrWDHQtgNAiQwMgw8gww2awdDKPPOpDAyDjfyDMAu5+1g2H0CZDAyDDyDDDZ5udm05vas67N2sEwWgRIYGQYeQaYbAf3z+TIoX2Zme6lksxM93Lk0D63McAIcQ8kMDLm52bX3QOZGHkGmDTWDobRJkACI2P1hEEVVgCA0SRAAiPFyDMAwOgSIAEA2FHW+IXdQ4AExpqTEoDRZo1f2F1UYQXG1upJydLySlruPilZWFwadtcAGLDGL+wuAiQwtpyUAIw+a/zC7iJAAmPLSQnA6LPGL+wuAiQwtpyUAIy++bnZ9Kb2rGuzxi+MLwESGFtOSgBG38H9MzlyaF9mpnupJDPTvRw5tE8BHRhTqrACY2v15EMVVoDRZo1f2D0ESGCsOSkBALh0BEgAADbFGrwwuQRIYGI44QG4eKtr8K4uo7S6Bm8Sn6kwARTRASbC6gnP0vJKWu4+4VlYXBp21wDGijV4YbIJkMBEcMIDsD2swQuTTYAEJoITHoDtYQ1emGwCJDARnPAAbA9r8MJkEyCBieCEB2B7HNw/kyOH9mVmupdKMjPdy5FD+xTQgQmhCiswEVZPbFRhBbh41uCFySVAAhPDCQ/AxixzBHQlQAKcgxMqYBJY1xHYDPdAAmzAupHApLDMEbAZAiTABpxQAZPCMkfAZgiQABtwQgVMCsscAZshQAJswAkVMCkscwRshgAJsAEnVMCksK4jsBmqsAJswLqRwCSxzBHQlQAJcA5OqIBxZRkiYKcIkADbyEkbMGzWdQR2knsgAbaJtSOBUWAZImAnCZAA28RJGzAKLEME7CQBEmCbOGkDRoFliICdJEACbBMnbcAosAwRsJMESIBt4qQNGAXWdQR2kiqsANvE2pHATthKdWfLEAE7RYAE2EabPWmz7AdwPpbkAEaNKawAQ2LZD+BCVHcGRo0ACTAkTgyBC1HdGRg1AiTAkDgxBC5EdWdg1AiQAEPixBC4ENWdgVEjQAIMiRNDmEwLi0t5yg1vzcOvf0OecsNbz3vfsyU5gFGjCivAkGx12Q+VW2F8baWqqiU5gFEiQAIM0VaW/VDSH8bX+Ypn+TcMjANTWAHGiMqtMN4UzwLGnQAJMEacfMJ4UzwLGHdDCZBV9ayqen9V3VlVB9a0P6yqVqrqvYPHr6957vFVdaKqbqmqX66qGrTfp6reXFU3D75+2aC9BtvdUlV/WVWPu/R/KcD2cvIJ403xLGDcDesK5PuSHEry9g2e+1Br7asHj+etaf+1JM9Ncu3g8bRB+/VJ3tJauzbJWwY/J8m3rtn2uYPXA4y1rZ58bqbqI9DdZv9tqaoKjLuhFNFprX0gSQYXES+oqh6c5Etaa386+Pk3kxxM8sYkT0/yTwebvjLJHyf58UH7b7bWWpJ3VdV0VT24tXbb9v0lAJfWViq3KrwDO2Or/7ZUVQXG2ShWYX14VS0m+cck/7619j+TzCS5dc02tw7akuSBq6GwtXZbVT1g0D6T5KMbvOYeAbKqnpv+Vco89KEP3cY/BWD7bfbkU9VH2Bn+bQGTaMcCZFX9YZIHbfDUT7TWXneOl92W5KGttb+vqscnWaiqxyTZ6FJlu1AXur6mtfaSJC9JkgMHDlzofQHGisI7sDP82wIm0Y4FyNbaN2/hNV9I8oXB9++pqg8leVT6Vw8fsmbThyQ5Nfj+71anpg6mun580H5rki8/x2sAJsbe6V6WNjihVXgH1ltYXNrU9HD/toBJNFLLeFTV/atqz+D7r0y/AM5fD6aofrqqnjSovvq9SVavYr4+ybMH3z/7rPbvHVRjfVKST7n/EZhECu/Aha3ez7i0vJKWu+9nPN9+r6IqMImGtYzHM6rq1iRPTvKGqjo+eOrrk/xlVf1FktckeV5r7R8Gzz0/yUuT3JLkQ+kX0EmSG5I8tapuTvLUwc9JcmOSvx5s/1+T/Kud/asARtNWqj5u5WQaxtn57mc8FxVVgUlU/SKlrDpw4EC76aabht0NgKF6yg1v3XBq3sx0L++4/roh9Ah21sOvf8OGhRIqyYdv+PZL3R2Aoaqq97TWDmz03EhNYQVgNCgOwqQ5132L7mcEWG8Ul/EAYMi2Whxks0VIYKdsdl+cn5tdt6Zj4n5GgI24AgnAPWylOIj7JhkVW9kX3c8I0I0rkADcw+pJ82au4FhUnVGx1X3x4P4Z+yrABQiQAGxosyfT7ptkJ2xlWrR9EWDnCJAAbIuLWVTdvZNsZHUq6urVxNWpqEnOu39czL4IwPm5BxKAbbHVRdXdO8m5bGVtxmTr+yIAF+YKJADbYiv3Ta5u797JybDZK81bnYq61X0RgAsTIAHYNlspQrLVkGDa63jZynTUi5mKqiAOwM4whRWAodrKAu6mvY6frUxHNRUVYPQIkAAM1VZCwlbvjUv64fMpN7w1D7/+DXnKDW8VOrdgK/8Nt3Kl2dqMAKPHFFYAhmor96tdzLTXrVT13M02OxX4UldGNRUVYLQIkAAM3WZDwlbDyFYL9ozD/ZZb6eNWwuBW/xvOz82u+12J6agA40iABGDsbDWMbOXK5cVctdxqqLsUQTDZWhhUGRVgsgmQAIydrYaRrVy5vJirlpsNdZcyCCZbC4MqowJMNkV0ABhLB/fP5B3XX5cP3/Dtecf113UKJlsp2LPVK25bKfSz1eJAW+3jVirgqowKMNkESAAmxlaqem4lZCVbC3WXMggmWwuDKqMCTDZTWAGYKJudRrnV+y23MtVzq9NDt9rHrU4FNhUVYHIJkABwHlsNWVsJdZc6CK6+VhgEoKtqrQ27DyPlwIED7aabbhp2NwDYBS5VFVYA2E5V9Z7W2oENnxMg1xMgAQCASXa+AKmIDgAAAJ0IkAAAAHQiQAIAANCJAAkAAEAnAiQAAACdCJAAAAB0IkACAADQiQAJAABAJwIkAAAAnQiQAAAAdCJAAgAA0IkACQAAQCcCJAAAAJ0IkAAAAHQiQAIAANCJAAkAAEAnAiQAAACdCJAAAAB0IkACAADQiQAJAABAJwIkAAAAnQiQAAAAdCJAAgAA0IkACQAAQCcCJAAAAJ0IkAAAAHQiQAIAANCJAAkAAEAnAiQAAACdCJAAAAB0IkACAADQiQAJAABAJwIkAAAAnQiQAAAAdCJAAgAA0IkACQAAQCcCJAAAAJ0IkAAAAHQiQAIAANCJAAkAAEAnAiQAAACdCJAAAAB0IkACAADQiQAJAABAJwIkAAAAnQiQAAAAdCJAAgAA0IkACQAAQCcCJAAAAJ0IkAAAAHQiQAIAANCJAAkAAEAnAiQAAACdCJAAAAB0IkACAADQiQAJAABAJwIkAAAAnQiQAAAAdCJAAgAA0IkACQAAQCcCJAAAAJ0IkAAAAHQiQAIAANCJAAkAAEAnAiQAAACdCJAAAAB0IkACAADQiQAJAABAJwIkAAAAnQiQAAAAdCJAAgAA0IkACQAAQCcCJAAAAJ0IkAAAAHQiQAIAANCJAAkAAEAnAiQAAACdCJAAAAB0IkACAADQiQAJAABAJwIkAAAAnQiQAAAAdCJAAgAA0IkACQAAQCcCJAAAAJ0IkAAAAHQiQAIAANCJAAkAAEAnAiQAAACdCJAAAAB0IkACAADQiQAJAABAJ9VaG3YfRkpVfSLJ3wy7HwzN/ZLcPuxOMNbsQ2wH+xEXyz7ExbIPTbavaK3df6MnBEhYo6puaq1m9ggKAAAIrUlEQVQdGHY/GF/2IbaD/YiLZR/iYtmHOBdTWAEAAOhEgAQAAKATARLWe8mwO8DYsw+xHexHXCz7EBfLPsSG3AMJAABAJ65AAgAA0IkACQAAQCcCJLtaVT2rqt5fVXdW1YGznjtcVbdU1cmqmlvT/rRB2y1Vdf2a9odX1Z9V1c1V9TtVdcWg/V6Dn28ZPP+wS/X3celV1U9V1VJVvXfw+LY1z23LPsXkOte+AklSVR+pqhODz56bBm33qao3Dz5H3lxVXzZor6r65cG+9JdV9bg17/PswfY3V9Wzh/X3cGlU1cur6uNV9b41bdu231TV4wf75S2D19al/Qu51ARIdrv3JTmU5O1rG6vq0Um+K8ljkjwtya9W1Z6q2pPkvyT51iSPTvLdg22T5OeT/FJr7dokn0zynEH7c5J8srX2yCS/NNiO3e2XWmtfPXjcmGz7PsUEusC+Aqu+cfDZszooen2Stww+R94y+Dnp70fXDh7PTfJrST84JPkPSZ6Y5AlJ/sNqeGDXekX6x6W1tnO/+bXBtquvO/t3scsIkOxqrbUPtNZObvDU05O8qrX2hdbah5Pckv4H4hOS3NJa++vW2heTvCrJ0wejadclec3g9a9McnDNe71y8P1rknyT0beJtJ37FJNpw31lyH1i9K09Bp19bPrN1veuJNNV9eAkc0ne3Fr7h9baJ5O8OU74d7XW2tuT/MNZzduy3wye+5LW2p+2fmXO34xj2a4nQDKpZpJ8dM3Ptw7aztV+3yTLrbU7zmpf916D5z812J7d64WDqT0vXzMCu537FJPpXPsKrGpJ/qCq3lNVzx20PbC1dluSDL4+YNC+2c8kJst27Tczg+/PbmcXu3zYHYCLVVV/mORBGzz1E621153rZRu0tWw8qNLOs/353osxdb59Kv2pOi9K///xi5L8pyT/e7Z3n2Iy2Se4kKe01k5V1QOSvLmqPniebc+1P9nPOJ/N7jf2pwkkQDL2WmvfvIWX3Zrky9f8/JAkpwbfb9R+e/rTOC4fXDFau/3qe91aVZcn+dLcc6oIY6TrPlVV/zXJ7w9+3M59isl0vn0I0lo7Nfj68ap6bfrTnv+uqh7cWrttMJ3w44PNz7U/3Zrkn57V/sc73HVGz3btN7cOvj97e3YxU1iZVK9P8l2DCqoPT/+m7/+V5N1Jrh1Ux7wi/aIorx/M6/+jJM8cvP7ZSV635r1Wq5E9M8lbB9uzCw0OtKuekX6hpmR79ykm04b7ypD7xIioqqur6t6r3yf5lvQ/f9Yeg84+Nn3voKrmk5J8ajBV8XiSb6mqLxtMwf+WQRuTZVv2m8Fzn66qJw3u7f/eOJbteq5AsqtV1TOS/EqS+yd5Q1W9t7U211p7f1W9OslfJbkjyQtaa2cGr3lh+h+Ue5K8vLX2/sHb/XiSV1XVi5MsJnnZoP1lSf5bVd2S/pXH77pEfx7D8QtV9dXpT9H5SJIfTJJt3qeYQK21O86zr8ADk7x2UKPt8iS/3Vp7U1W9O8mrq+o5Sf42ybMG29+Y5NvSL+j1uSTfnySttX+oqhelP2CRJD/TWjNrZherqv+R/tXD+1XVrelXU70h27ffPD/9Sq+9JG8cPNjFyoUSAAAAujCFFQAAgE4ESAAAADoRIAEAAOhEgAQAAKATARIAAIBOBEgAuASq6nlV9b2D77+vqvauee6lVfXobfo9B6vqJwffv6Kqnnmh15zjfe5fVW/ajj4BsHtYBxIALoHW2q+v+fH70l8E/tTguf9jG3/VjyX5Zxf7Jq21T1TVbVX1lNbaO7ahXwDsAq5AAjCxquprquovq+rKqrq6qt5fVV911jYPq6oPVtUrB9u+pqquGjz3TVW1WFUnqurlVXWvQfsNVfVXg+3/46Dtp6rq/xxcETyQ5Leq6r1V1auqP66qA4Ptvnvwfu+rqp9f04/PVNXPVtVfVNW7quqBG/w9j0ryhdba7Rs896LBFcnLquojVfVzVfWnVXVTVT2uqo5X1Yeq6nlrXraQ5Hsu+j80ALuGAAnAxGqtvTvJ65O8OMkvJPnvrbX3bbDpbJKXtNb+tyT/mORfVdWVSV6R5Dtba/vSn9Xz/Kq6T5JnJHnMYPsXn/U7X5PkpiTf01r76tbayupzg2mtP5/kuiRfneRrqurg4Omrk7yrtfbYJG9P8gMb9PMpSf787Maq+oUkD0jy/a21OwfNH22tPTnJ/xz8Hc9M8qQkP7PmpTcl+boNfg8AE0qABGDS/UySp6Z/VfAXzrHNR9dM4/zvSb42/VD54dba/zdof2WSr08/YH4+yUur6lCSz22iL1+T5I9ba59ord2R5LcG75kkX0zy+4Pv35PkYRu8/sFJPnFW2/+VZLq19oOttbam/fWDryeS/Flr7dOttU8k+XxVTQ+e+3iSvQGAAQESgEl3nyTXJLl3kivPsU3b4OfacMN+8HtCkt9LcjDJZgrRbPieA6fXBMAz2biOwUru+Te8O8njB1dG1/rC4Ouda75f/Xn1va8cvCcAJBEgAeAl6V+l+630p49u5KFV9eTB99+d5E+SfDDJw6rqkYP2f5nkbVV1TZIvba3dmOSH05+KerZPpx9Yz/ZnSb6hqu5XVXsGv+ttm/hbPpDkkWe1vSnJDUneUFUb/c7zeVT6xX4AIIkqrABMsMGyGne01n57ENjeWVXXtdbeetamH0jy7Kr6jSQ3J/m11trnq+r7k/xuVV2e/pW+X0//iubrBvdIVpIf2eBXvyLJr1fVSpLVYJrW2m1VdTjJHw1ee2Nr7XWb+JPenuQ/VVWtna7aWvvdQXh8fVV92ybe7xuTvGET2wOwy9X62yEAgLWq6mFJfr+19lUX2HQkVNX/neT/ba394Ta819uTPL219smL7xkAu4EprACwu/xckqsu9k2q6v5JflF4BGAtVyABAADoxBVIAAAAOhEgAQAA6ESABAAAoBMBEgAAgE4ESAAAADr5/wF/FTODZUtlhwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x1080 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(15, 15))\n",
    "\n",
    "x_positions = positions[:, 0]\n",
    "y_positions = positions[:, 1]\n",
    "plt.scatter(x_positions, y_positions)\n",
    "\n",
    "plt.xlabel(\"x position (km)\")\n",
    "plt.ylabel(\"y position (km)\")\n",
    "\n",
    "# Some gibberish that draws a big blue circle representing the Earth:\n",
    "earth = plt.Circle((0, 0), radius_of_earth, color='b')\n",
    "plt.gcf().gca().add_artist(earth)\n",
    "\n",
    "# Make the plot big enough to show elliptical orbits:\n",
    "plot_limit = 8000\n",
    "plt.xlim(-1.8 * plot_limit, 1.8 * plot_limit)\n",
    "plt.ylim(-2.4 * plot_limit, 1.2 * plot_limit)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Playing with the Notebook and Checkout Questions\n",
    "\n",
    "### Question 1a\n",
    "\n",
    "Using the intial velocity of 8.982 km/s, where in the ellipse is the cannonball moving the slowest? How are you answering this?\n",
    "\n",
    "### Question 1b\n",
    "\n",
    "Where in the ellipse is the cannonball moving the fastest?\n",
    "\n",
    "### Question 2a\n",
    "\n",
    "Without changing the simulation duration, increase the initial velocity. Do you expect to get a larger or smaller ellipse? Is the ellipse more or less eccentric (more oval or more circular)?\n",
    "\n",
    "### Question 2b\n",
    "\n",
    "With the increased initial velocity, do you expect the orbit to over-close in the same simulation duration, or will it not close?\n",
    "\n",
    "### Question 3a\n",
    "\n",
    "Still without changing the simulation duration, decrease the initial velocity. Do you expect to get a larger or smaller ellipse? Is the ellipse more or less eccentric (more oval or more circular)?\n",
    "\n",
    "### Question 3b\n",
    "\n",
    "With the decreased initial velocity, do you expect the orbit to over-close in the same simulation duration, or will it not close?\n",
    "\n",
    "### Question 4\n",
    "\n",
    "What value of the initial velocity gets a circle?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "With oodles of software and hardware and just three weeks of time learning to use it, you can get more results than Isaac Newton."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}