Classical_mechanics_solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Instructions\n",
    "\n",
    "Read carefully all the content in the cells, including the comments in the code cells\n",
    "\n",
    "Some lines of code have been removed. You are asked to fill in the gaps.\n",
    "\n",
    "Only modify lines that are between the following promts:\n",
    "\n",
    "'# START YOUR CODE HERE >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n",
    "\n",
    "'# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> END YOUR CODE HERE\n",
    "\n",
    "After each cell there will be a test cell that will allow you to be confident you will be able to continue with the execution of the notebook.\n",
    "\n",
    "Here's an example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:37:15.565541Z",
     "start_time": "2021-01-31T23:37:15.561326Z"
    }
   },
   "outputs": [],
   "source": [
    "# This is just a test to illustrate how we will do exercises with jupyter notebooks\n",
    "\n",
    "# This function adds two numbers and prints the result\n",
    "# Please fill in the gaps\n",
    "\n",
    "def print_sum(a, b):\n",
    "\n",
    "    # START YOUR CODE HERE (1 line of code)>>>>>>>>>>>>>>>>>>\n",
    "    c = a + b\n",
    "    # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> END YOUR CODE HERE\n",
    "\n",
    "    print(c)\n",
    "\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:37:17.111453Z",
     "start_time": "2021-01-31T23:37:17.106661Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3\n"
     ]
    }
   ],
   "source": [
    "# Run this cell to check that the above implementation is correct\n",
    "print_sum(1, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Expected Output**: \n",
    "<table style = \"width:40%\">\n",
    "    <tr>\n",
    "    <td>** print_sum(1,2) **</td> \n",
    "        <td>3 </td> \n",
    "    </tr>\n",
    "\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Newton's laws of motions\n",
    "\n",
    "In this notebook, we are going to keep practicing how to code in Python, and use notebooks as a coding / development environment, to solve simple problems using Newton's laws of motion.\n",
    "\n",
    "In particular, we are going to see how to solve a problem you may have already encountered in school: figuring out which angle to throw a ball at to maximize the distance it reaches when falling on the ground.\n",
    "\n",
    "These are the simplifying hypothesis we will use:\n",
    "\n",
    "1. The ball's dimensions and geometry are not relevant for the problem (it is very small compared to the distance travelled): we will consider it to be a point particle with mass $m$.\n",
    "\n",
    "2. The Earth is flat (we are not launching a rocket with a range comparable to its radius). We will also consider that the force of gravity is constant, with its acceleration equal to $g=9.81\\,\\rm m \\, s^{-2}$.\n",
    "\n",
    "3. There is no atmosphere, so we will not consider the effect of air resistance.\n",
    "\n",
    "4. The ball is initially launched with a velocity $v_0$ at an angle $\\theta$ from the horizontal.\n",
    "\n",
    "Under these assumptions, the motion of the ball is determined by Newton's second law of motion:\n",
    "\n",
    "$$\\mathbf{F} = m \\mathbf{a}$$\n",
    "\n",
    "Let's set up a coordinate system, with the $x$ axis parallel to the ground and pointing towards the initial velocity of the ball, and the $y$ axis perpendicular to the ground pointing upwards. The $z$ axis is set so that $\\{x,y,z\\}$ is a right-handed coordinate system. The components of the equations of motion are then:\n",
    "\n",
    "$$a_x = 0 \\\\ a_y = -g \\\\ a_z = 0$$\n",
    "\n",
    "If we solve these equations, we can find the ball's coordinates as a function of time:\n",
    "\n",
    "$$\\{x(t), y(t), z(t)\\}$$\n",
    "\n",
    "From which we can compute the distance traveled by the time the ball hits the ground. \n",
    "\n",
    "We will follow two approaches:\n",
    "\n",
    "  a. This system of equations has an analytic solution. We will compute it and use it to solve the problem.\n",
    "  b. We will solve it numerically, following an approach that can be used for more general problems, for which there is no analytic solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## a. Analytic solution\n",
    "\n",
    "Let's use the definition of acceleration and velocity to integrate the equations of motion. First, integrating the acceleration we get the velocity (since the acceleration is the derivative of the velocity):\n",
    "\n",
    "$$a_x \\equiv \\frac{d v_x}{d t} = 0 \\\\ a_y \\equiv \\frac{d v_y}{d t} = -g \\\\ a_z \\equiv \\frac{d v_z}{d t} = 0 $$\n",
    "\n",
    "The derivative of a straight line is its slope:\n",
    "\n",
    "$$\\frac{d}{d x} \\left(a x + b \\right) = a$$\n",
    "\n",
    "Then:\n",
    "\n",
    "$$ v_x = v_{x0} \\\\ v_y = -g t + v_{y0} \\\\ v_z = v_{z0} $$\n",
    "\n",
    "Since there is no force along the $y$ and $z$ axes (due to our choice of reference frame), there is no acceleration, and the velocity of the ball does not change (it keeps the same value it was launched with). A constant force (in this case, along the $y$ axis) induces a linear change in te velocity.\n",
    "\n",
    "Now, we can integrate again to find the position:\n",
    "\n",
    "$$ v_x \\equiv \\frac{d x}{d t} \\rightarrow x = v_{x0} t + x_0 \\\\\n",
    "   v_y \\equiv \\frac{d y}{d t} \\rightarrow y = -\\frac{g}{2} t^2 + v_{y0} t + y_0 \\\\\n",
    "   v_z \\equiv \\frac{d z}{d t} \\rightarrow z = v_{z0} t + z_0$$\n",
    "\n",
    "Constant velocities induce linear changes in positions. Linear velocities induce quadratic changes in positions. \n",
    "\n",
    "The initial velocity of the ball is:\n",
    "\n",
    "$$v_{x0} = v_0 \\cos \\theta \\\\ v_{y0} = v_0 \\sin \\theta \\\\ v_{z0} = 0$$\n",
    "\n",
    "The initial position of the ball is:\n",
    "\n",
    "$$ \\{x_0, y_0, z_0\\} = \\{0,0,0\\}$$\n",
    "\n",
    "Substituting, the equations of motion are:\n",
    "\n",
    "$$x = v_0 \\cos (\\theta) t \\\\ y = -\\frac{g}{2}t^2 + v_0 \\sin (\\theta) t \\\\ z = 0$$\n",
    "\n",
    "We can compute the position of the ball with time, and plot it using Python.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:37:21.984860Z",
     "start_time": "2021-01-31T23:37:21.201154Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAFgCAYAAABEyiulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3iV9f3/8ec7ISHMBEgggbDDFkGICIKIrYuhVjvUn4p1IVX71dY66qjWWrXLqrVqVVBRFLdFZbgVVJA9EyDsEDJYGZCdz++PHGxEwAAn5z7j9biuXJxznzv3eR24yCufe3xuc84hIiIioSXK6wAiIiJy5FTgIiIiIUgFLiIiEoJU4CIiIiFIBS4iIhKCGnkdwB8SExNdly5dvI4hIiLid4sWLdrhnEs6cHlYFHiXLl1YuHCh1zFERET8zsw2H2y5dqGLiIiEIBW4iIhICFKBi4iIhCAVuIiISAhSgYuIiIQgFbiIiEgIUoGLiIiEIBW4iIhICFKBi4iIhCAVuIiISAhSgYuIiISgsJgLXUQkmDjnKK2spryyhorqGiqq6vx5wOPyA1+rqv7u63UeR0cZ7VrG0T4hjpT4JrSPb0K7+MY0bhTt9UcWD6jARUSOknOOvKJy1uYVszavmKz8EtbmFbMuv4Tisiq/vEdsdBSxjWq/KqtrDrrdxOaxpMQ3ITk+jvbxcaQkNCElvrbkU+LjSI6PIyZaO1zDjQpcROQHOOfILSpjXV7JYYu6dbNYerRtznkD25PaqimNfcW7v4QbN4oiJvq7y/Yvj42O/vb5/tdjog0z+06WveVVbC8sY3thae2fe/73ePPOvcxbv5Pi8u+WvBkkNm9cW+7xTUhJiPu24Du1bkrf9i1V8CEooAVuZh2BKUAyUAM87Zx79IB1RgH/BTb6Fr3lnLsvkDlFJDLtL+q1eSWsyytmXV4J6/IPXdQ/GdiBHu2a06NtC3q2a06b5o0bPGOzxo1Ia9uctLbND7lOcVkluYVl5BSWkVtYSk6dks8qKGHOugL2VlT/b5ux0aR3ac3Qbm0Y2q01/TvE00iFHvQCPQKvAm52zi02sxbAIjP70Dm3+oD15jjnxgU4m4hEmILicj7NzGfR5t2szS8mK6/kO6PXNs1iSfMVdc92zUkLYFEfixZxMbSIi6FHuxYHfd05R3F5Fdv3lLEuv5j5G3Yxb8NO/jIrE6gt9BO77i/0NhzXvqUKPQgFtMCdc9uB7b7HxWaWAXQADixwERG/c86RlV/Chxl5fLQ6jyVb9+BcnRH1CaFV1EfLzGgZF0PL5Bh6Jbdg3PHtgdpfaOZv3Mm8DTuZt2EXD82sLfTmjRtxYpdWDO3WhmHd29A3RYUeDMw5580bm3UBvgCOc84V1Vk+CngTyAZygN8551Yd5PsnABMAOnXqNHjz5s0NH1pEQk5VdQ0LN+/mo9V5fJiRx+ad+wDo3yGeM/q24/Q+7eiT0uJ7x5oF8ovLvh2df71hJxsK9gLQonEj3wi9NcO6JdK3fUuio/T311DMbJFzLv17y70ocDNrDnwO/Nk599YBr7UEapxzJWY2BnjUOdfjcNtLT093CxcubLjAIhJSSsqr+HxNAR9l5PFJZj6FpZXERkdxclobTu/Tjh/3aUtKfBOvY4ac/KIy5m3cxdfrdzJ/w0427PhfoQ/p2pph3Wt3ufdJUaH7U9AUuJnFAO8Bs51zD9dj/U1AunNux6HWUYGLSM6eUj7OyOPDjHzmrd9JRXUNCU1j+FHvtpzRpx2n9EyieWNdeONPeUVlvt3ttbvcN/oKvVXTGM4b2IGLhnSkd3JLj1OGvqAocKvdR/UCsMs5d9Mh1kkG8pxzzsyGAG8And1hgqrARSKPc45VOUV8lJHHh6vzWJVTeySuS5umnNG3HWf0TWZQpwQdqw2g3MLaQv8wI48PV+VRUV3DgI4JXHRiR84Z0F6/QB2lYCnwEcAcYAW1l5EB3AF0AnDOPWVmNwC/ovaM9VLgt865rw63XRW4SGSoqKph3oadfOQ7CS2nsAwzGNypFaf7jmd3T2qm49lBYNfeCt5eso1XF2xhbV4JTWOjGXd8Chee2IlBnRL0b3QEgqLAG4oKXCS87Swp56V5W3hx3iZ2lFTQJCaaU3okcnrfdvyod1sSw/Rs8XDgnGPJ1j28+s1W3l2ew76Kanq0bc6FJ3bkgkGptG4W63XEoKcCF5GQsy6vmMlfbuTNxduoqKrhtF5JXHJSZ0b0SCQuRvN/h5qS8ireX57DtAVbWbJlDzHRxpn9krnoxI4M755IlE58OygVuIiEBOccc7N2MGnuRj5bU0DjRlFcMCiVq0Z0Ia3twScmkdCzJreYVxds5a0l2ezZV0lqqyb8Ir0jP09P1RUCB1CBi0hQK6+q5r9Lc5g8dyOZucUkNm/M+GGdueSkTmE7oYrU/rt/sCqPVxdsZW7WDqIMTu2ZxIUnduLHfdpqjnZU4CISpHaWlDN1/hamfL2ZHSXl9E5uwVUjunLuwPa6TWaE2bprH68t3MrrC7PJLSojsXksPx2cyoXpHemWdOi538OdClxEgkpWfjGT5m7ircXZlFfVMKpXEleP6MbwtDY6QznCVdc4vlhbwLQFW/g4I5+qGseQrq258cc9GJ6W6HW8gFOBi4jnnHN8mbWTZ+du0PFtqZf84jLeWryNKV9tIqewjOFpbbj1rN4M6JjgdbSAUYGLiGfKq6qZvjSHSd8e345l/LAuOr4t9VZeVc3UeVt4/NMsdu2tYPRxydx8Zq/D3lY1XKjARSTgdu2tYOq8zbzgO77dq10LrjqlK+cOaK/LwOSolJRXMWnORp6Zs4F9FVX8bHAqN57ekw4J4XvmugpcRAJm994KHv14Ha98s+Xb49tXjejKiLREHd8Wv9hZUs4Tn63nxa9r70R52bDOXDeqe1ju0VGBi0iDq6qu4ZVvtvCPD9dSXFbFzwalcvUpXenRTse3pWFs21PKox+t5Y1F2TSJieaakd24+pRuYTXvugpcRBrU1+t38sd3V5GZW8zJ3dtwzzn96JWs4pbAyMov5h8frGXmylxaN4vl+tPSuOSkTmFxqEYFLiINInv3Ph6ckcn7K7bTIaEJd43tw9nHJWtXuXhi2dY9/G32GuZm7aBDQhNuPL0HF5zQIaTvSqcCFxG/Kqus5qnP1/PkZ+sxg+tGpTFhZLewGPFI6Psyawd/nZXJsuxC0to253dn9uSsfqH5i6UKXET8wjnHzJW5/Pn9DLbtKWXc8Sn8fkyfsD4LWEKTc47Zq3L52+w1rC/Yy4COCdx6Vq+QmwxGBS4ixywzt4h7p69i3oZd9E5uwb3n9mNotzZexxI5rKrqGt5aso1HPlxLTmEZI9ISueWsXiEzGYwKXESO2p59FTz84VpemreZlk1iuPnMXlx8YseQPq4okaesspqp87fw7zqTwdw1rm/Q7z1SgYvIEauucbz8zRb+8cEaikoruXRoZ357Rk8SmsZ6HU3kqBWXVTJp7kae/mIDUWbcNbYPF57YMWiPj6vAReSIzN+wk3vfXU3G9iKGdmvNPef0o09KS69jifjN1l37uOWNZczbsIuRPZN46IL+tA/C0bgKXETqJWdPKQ/MyOC95bWXhd05tg+jdVmYhKmaGsdL8zfz4IxMGkUZd43rwy/Sg2s0rgIXkcMqq6zm6S828MRnWTgHvxrVnWtHdqdJrC4Lk/C3ZWftaHz+xl2c2jOJh37an5T44BiNq8BF5JA+WJXLfe+tJnt3KWP7p/D7Mb1JbdXU61giAVVT43hx3mYemplJo2jj7nF9+fngVM9H44cq8PCZLFZEjti+iirue3c10xZspXdyC165ZijDuuuyMIlMUVHG5Sd3YVSvJG55fTm3vrGcmSu28+AFx5McH+d1vO/RCFwkQq3KKeT/XlnChh17+dWp3fnNGT2J0WVhIkDtaPyFrzfxl1mZxERHcc85/fjpoA6ejMYPNQLX/1aRCOOcY9LcjZz/768oKa9i6lUncevZvVXeInVERRlXDO/KrBtH0ju5Bb97fRlXvbCQvKIyr6N9S/9jRSLIjpJyrnx+AX96bzUjeyYx88aRnBxi00qKBFKXxGa8OmEYd4/ry1frd3DGw5/z1uJsgmHvtQpcJEJ8sbaAsx+Zw5frd/Kn8/rxzPjBtG6mCVlEfkhUlHHViK7MvHEkPdu14LevLeOaKQvJ93g0rgIXCXMVVTU8MCOD8ZO/oXWzGKbfMJzLhnXx/MxakVDTNbEZr147jLvG9mHOuh2c8c8veGfJNs9G4ypwkTC2ccdefvrkVzz9xQYuHdqJ6TeMoHeyZlMTOVrRUcbVp3Rjxo2n0D2pGTe9upQJLy4ivzjwo3EVuEgYcs7x+sKtjH1sDlt37+M/lw3m/p/01726Rfyke1JzXp94MneO6cMXaws4859f8N+lgR2Nq8BFwkxRWSX/N20pt7yxnONT45l54ymc1S/Z61giYSc6yrhmZO1ovGtiM26ctpRrX1xEQXF5QN5fBS4SRhZt3s2YR+cwY8V2bjmrF1OvHho000GKhKvuSc15Y+LJ3DGmN5+tLeCbjbsC8r6aiU0kDFTXOJ74NItHPl5H+4Q4Xp84jEGdWnkdSyRiREcZE0Z255wB7QP2S7MKXCTEbS8s5aZpS5m/cRfnDmjP/ecfR8u4GK9jiUSkQO7xUoGLhLBZK3O57c3lVFbX8PefD/BsqkcRCTwVuEgIKq2o5v73VzN1/hb6d4jnsYtPoGtiM69jiUgAqcBFQkxmbhG/fnkJ6/JLmDCyG787sxexjXQ+qkikUYGLhJBZK3O56dUlNG8cw5QrhzCyZ5LXkUTEIypwkRCw/w5if56RwYDUBJ4Zn05Si8ZexxIRD6nARYJcVXUNf3x3NS/O28zo45L554UDNaOaiKjARYJZSXkVv355MZ+uKeDakd247ezeREXpLHMRUYGLBK3cwjKufH4Ba/KK+fP5x3HJSZ29jiQiQUQFLhKEVucUceXzCyguq2TS5emM6tXW60giEmRU4CJB5tM1+dwwdTEt4mJ4feLJ9G2v23+KyPepwEWCyEvzNnPP9FX0Tm7BpMtPJDk+zutIIhKkVOAiQaCmxvHgzAyembORH/Vuy78uPoFmjfXfU0QOLaDTN5lZRzP71MwyzGyVmd14kHXMzB4zsywzW25mgwKZUSTQSiuquW7qYp6Zs5Hxwzrz9GWDVd4i8oMC/VOiCrjZObfYzFoAi8zsQ+fc6jrrjAZ6+L5OAp70/SkSdgqKy7l6ykKWZ+/h7nF9uXJ4F92MRETqJaAF7pzbDmz3PS42swygA1C3wM8DpjjnHDDPzBLMLMX3vSJhIyu/mF8+t4AdJeU8delgzuqX7HUkEQkhnu2nM7MuwAnA/ANe6gBsrfM827fsOwVuZhOACQCdOnVqqJgiDeKr9TuY+OIiYhtF89q1wzg+NcHrSCISYjy5hZGZNQfeBG5yzhUd+PJBvsV9b4FzTzvn0p1z6UlJuqGDhI43FmUzftI3JMfH8c71J6u8ReSoBHwEbmYx1Jb3VOfcWwdZJRvoWOd5KpATiGwiDck5xz8/XMtjn2QxIi2RJy4dRMu4GK9jiUiICvRZ6AZMAjKccw8fYrXpwHjf2ehDgUId/5ZQV15VzU2vLuWxT7K4ML0jz11xospbRI5JoEfgw4HLgBVmttS37A6gE4Bz7ilgBjAGyAL2AVcEOKOIX+3eW8G1Ly7im027uOWsXlw3qrvONBeRYxbos9DncvBj3HXXccD1gUkk0rA27djLFc8vYNueUv518QmcM6C915FEJExotgiRBpKxvYhLnp2Pc46Xrz6J9C6tvY4kImFEBS7SAPaXd+NGUbx8zVC6JjbzOpKIhBkVuIifZeb+r7xfuWYoXVTeItIAPLkOXCRcZeYW8f+emU9stMpbRBqWClzET+qW97QJKm8RaVgqcBE/UHmLSKCpwEWO0ZrcYpW3iAScClzkGKzJLebiZ+bVHvNWeYtIAKnARY7SgeWtS8VEJJBU4CJHoXa3+Txiok3lLSKeUIGLHKH95d0o2pg2YZjKW0Q8oQIXOQJr81TeIhIcVOAi9bQ2r5iLn64t71c0PaqIeEwFLlIPB5Z3t6TmXkcSkQinAhf5Aft3m0dHqbxFJHiowEUOY395R5kxbYLKW0SChwpc5BDWqbxFJIipwEUOYl1e7SQtKm8RCVYqcJED1C3vV1TeIhKkVOAiddSW9/xvy7u7yltEgpQKXMRnf3mbofIWkaCnAhcBNu/c+215T1N5i0gIUIFLxNu1t4JfPreAqpoaXrlG5S0ioaGR1wFEvFRWWc01UxaybU8pL199EmltVd4iEho0ApeIVVPj+O1rS1m8ZTePXjiQ9C6tvY4kIlJvKnCJWA/MyGDGilzuHNOH0f1TvI4jInJEVOASkZ77ciPPzt3IFcO7cPUp3byOIyJyxFTgEnFmrczlvvdWc1a/dtw1tq/XcUREjooKXCLKos27uXHaEgZ2TODRi04gOsq8jiQiclRU4BIxNu7Yy9UvLCAlPo5nx6cTFxPtdSQRkaOmApeIsLOknCue+wYz4/krhtCmeWOvI4mIHBMVuIS90opqrp6ykO2FZTx7eTpdEpt5HUlE5JhpIhcJa9U1jpteXcLSrXt48pLBDOrUyutIIiJ+oRG4hC3nHH96bzWzV+Xxh3F9Ofu4ZK8jiYj4jQpcwtakuRt5/qtNXD2iK1cM7+p1HBERv1KBS1iasWI7f56RwZj+ydwxpo/XcURE/E4FLmFn4aZd3PTqUgZ3asXDvxhIlK71FpEwpAKXsLKhoISrpywkNaEJz+habxEJYypwCRs7Ssr55XMLiPZd692qWazXkUREGowuI5OwsK+iiqueX0B+cRnTJgyjU5umXkcSEWlQGoFLyKuucfzfK0tZsa2Qf108iIEdE7yOJCLS4DQCl5DmnOPe6av4KCOP+87rxxl923kdSUQkIDQCl5D2zJwNvDhvM9eO7Mb4YV28jiMiEjAqcAlZ7y7L4YEZmYw7PoXbzu7tdRwRkYBSgUtI+mbjLm5+bRlDurTm7z8foGu9RSTiBLTAzWyymeWb2cpDvD7KzArNbKnv6w+BzCehYdOOvVwzZSEdWzfh6fGDda23iESkQJ/E9jzwODDlMOvMcc6NC0wcCTWlFdVMfGkRZvD8FUNIaKprvUUkMgV0BO6c+wLYFcj3lPDhnOOOt1ewJq+Yxy46gY6tda23iESuYDwGPszMlpnZTDPrd6iVzGyCmS00s4UFBQWBzCceeXHeZt5eso2bz+jJyJ5JXscREfFUsBX4YqCzc24A8C/gnUOt6Jx72jmX7pxLT0rSD/Nwt2jzLu57dzWn92nLdaPSvI4jIuK5oCpw51yRc67E93gGEGNmiR7HEo/lF5dx3dTFdGjVhH/o7mIiIkCQFbiZJZuZ+R4PoTbfTm9TiZcqq2u44eUlFJZW8tSlg4lvEuN1JBGRoBDQs9DN7BVgFJBoZtnAPUAMgHPuKeBnwK/MrAooBS5yzrlAZpTg8tdZmXyzcRf/vHAAfVJaeh1HRCRoBLTAnXMX/8Drj1N7mZkI7y/fzjNzNnL5sM6cf0Kq13FERILKIQvczA53rfbh3O2c23yU3ysCQFZ+Mbe8sYxBnRK4c2xfr+OIiASdw43ALwVygfIj2F4n4BFABS5HrbiskgkvLqJpbDRPXDKY2EZBdaqGiEhQ+KFd6D9xzn1Tnw2ZWSOg4tgjSSRzznHrG8vZvHMfL111EsnxcV5HEhEJSocb2nwMFB7Btqp931N0TIkkoj0zZwMzV+Zy+9m9Gda9jddxRESC1iFH4M65M45kQ76zxY/oe0Tq+mr9Dh6amcmY/slcfUpXr+OIiAQ1HVyUoLC9sJRfv7yEronN+OvPBuCbDkBERA6h3peR+Y5xDwM6At87MOmcm+zHXBJBKqpquG7qYsoqq/nPZUNp3jjQN8kTEQk99fpJaWaDgLeBVOBgQyMHqMDlqNz//mqWbNnDE5cMIq1tC6/jiIiEhPoOdZ4CSoCfAJnobHPxk7cWZzPl681MGNmNMf1TvI4jIhIy6lvgfYFf+G4wIuIXq3OKuOPtFQzt1ppbz+rldRwRkZBS35PY1gLNGjKIRJbCfZVMfGkR8U1i+NfFg2gUrfMpRUSORH1/at4B3GVmnRoyjESGmhrHb15byvbCUp64ZDBJLRp7HUlEJOTUaxe6c26WmY0C1pnZWmD391dxp/o7nISnxz/N4pPMfO47rx+DO7fyOo6ISEiq71notwO3AgXUzrRW3ZChJHx9tiaff360lvNP6MBlQzt7HUdEJGTV9yS2m4D/ADc451TeclS27trHjdOW0qtdCx44v78maxEROQb1PQbeFHhd5S1Hq6yyml9NXUSNczx16WCaxEZ7HUlEJKTVt8BnUjsLm8gRc85x9zsrWbmtiH/+YiBdEnVBg4jIsarvLvRHgOd9uzxn8f2T2HDObfBjLgkj0xZs5fVF2fz6R2mc3red13FERMJCfQv8S9+ffwLuO8Q62icq37NyWyH3/HcVp/RI5KbTe3odR0QkbNS3wK+kdr5zkXorrajmxmlLaNUshkcvOoHoKJ20JiLiL/W9Dvz5Bs4hYejBmRmsL9jLi1cNoXWzWK/jiIiEFc1fKQ3i08x8pny9matGdOWUHklexxERCTuHLHAzm2xmXeu7Ias1WdOtyo6Scm55Yzm9k1twi25SIiLSIA43Ar8cSDzCbR3p90iYcc5x+5srKCqt5JGLBhIXo3MbRUQawuGOgRvwpJkV1XNbOkNJmLZgKx9l5HHX2D70Tm7pdRwRkbB1uAL/gtozz4+kmL8Aio8pkYSsDQUl3PfuakakJXLl8HoffRERkaNwyAJ3zo0KYA4JcZXVNfzm1aXENori7z8fQJQuGRMRaVD1vQ5c5LAe+3gdy7ILeeKSQSTHx3kdR0Qk7OkyMjlmCzft4t+fZvGzwamM6Z/idRwRkYigApdjUlxWyW9eW0qHVk2455y+XscREYkY2oUux+Te6avZtruU1ycOo0VcjNdxREQihkbgctTeX76dNxdnc8NpaQzu3NrrOCIiEeVwM7H9OpBBJLTkFpZxx9srGNAxgV//uIfXcUREIs7hRuCPmNkXZpYWsDQSEmpqHDe/vpSKqhoeuXAgMdHakSMiEmiH+8n7YyAFWGZmvzMzXdgrAEz+ciNfZu3kD+f0pWtiM6/jiIhEpEMWuHPuM6A/8DjwAPC1mek04wiXsb2Iv85awxl923HRiR29jiMiErEOexa6c64MuM3MXgWeBRab2TtAxfdXdZc3UEYJEmWV1dw0bSktm8Tw0AX90U4ZERHv1PcysixgKTAQOIWDFLg/Q0lw+tvsNazJK+a5K06kTfPGXscREYloP1jgZnYu8ATQFLjGOTepwVNJ0JmzroBJczdy+bDOnNarrddxREQi3uEuI0vy7Tp/G1gMHKfyjky791bwu9eXkda2Ob8f08frOCIiwuFH4JlADTDeOTc1QHkkyDjnuOPtFezaW8Gky08kLiba60giIsLhLyP7BOin8o5sbyzKZubKXG4+sxfHdYj3Oo6IiPgc7n7gPw9kEAk+W3bu497pqxjarTXXnNLN6zgiIlKHptCSg6qqruE3ry0lKsr4xy8GEh2lS8ZERIKJ7kYmB/XEZ+tZtHk3j140kA4JTbyOIyIiBwjoCNzMJptZvpmtPMTrZmaPmVmWmS03s0GBzCe1lm7dw6Mfr+MnA9tz3sAOXscREZGDCPQu9OeBsw/z+migh+9rAvBkADJJHXvLq7hp2hKSW8bxx/OO8zqOiIgcQkAL3Dn3BbDrMKucB0xxteYBCWaWEph0AnD/+6vZvGsfD/9iAPFNYryOIyIihxBsJ7F1ALbWeZ7tW/Y9ZjbBzBaa2cKCgoKAhAt3n2Tm8co3W5l4andO6tbG6zgiInIYwVbgBzvV+aDzrDvnnnbOpTvn0pOSkho4VvgrLqvkzrdX0qtdC35zek+v44iIyA8ItrPQs4G696hMBXI8yhJR/jZ7DblFZTx56WBiGwXb73UiInKgYPtJPR0Y7zsbfShQ6Jzb7nWocLdw0y5enLeZK07uysCOCV7HERGRegjoCNzMXgFGAYlmlg3cA8QAOOeeAmYAY6i9fek+4IpA5otE5VXV3PbmctrHN+HmM7XrXEQkVAS0wJ1zF//A6w64PkBxBPj3p+tZX7CXF64cQrPGwXZERUREDiXYdqFLAK3JLebJz7K44IQOnNpTJwKKiIQSFXiEqq5x3PbmclrExXDXuL5exxERkSOkAo9QU77exNKte7jnnL60bhbrdRwRETlCKvAIlL17H3+bvYbTeiVx7oD2XscREZGjoAKPMM457ny79l4y95/fHzPdJlREJBSpwCPMf5fm8PnaAm49q5duEyoiEsJU4BFkZ0k5f3x3FSd0SuCyYV28jiMiIsdABR5B/vTeakrKq/jLT48nOkq7zkVEQpkKPEJ8uiafd5bmcN2oNHq2a+F1HBEROUYq8Aiwt7yKu95eSVrb5lx3Wnev44iIiB9o7swI8LfZa8gpLOWNiSfTuFG013FERMQPNAIPc4u37OaFrzcxfmhnBndu5XUcERHxExV4GKuoquH2N5eT3DKOW87u7XUcERHxI+1CD2NPfb6etXklTP5lOs11pzERkbCiEXiYysov5vFPsjh3QHt+1Lud13FERMTPVOBhqKbGcdubK2jaOJo/nKM7jYmIhCMVeBiaOn8zizbv5u6xfUls3tjrOCIi0gBU4GEmZ08pD83M5JQeiVwwqIPXcUREpIGowMOIc46731lJjYMHdKcxEZGwpgIPI+8t387HmfncfGZPOrZu6nUcERFpQCrwMLF7bwX3Tl/FgNR4rhje1es4IiLSwHRxcJi4//0MCksreenqk3SnMRGRCKAReBiYs66ANxdnM/HU7vRJael1HBERCQAVeIjbV1HFHW+voFtSM274UZrXcUREJEC0Cz3EPfzBWrbuKuW1a4cRF6M7jYmIRAqNwEPYsq17mPzlRi45qRNDurb2Oo6IiASQCjxEVVbXcNuby0lq0ZjbRutOYyIikUa70EPU5Lkbycwt5unLBtMyLsbrOCIiEiUZXg0AAA8VSURBVGAagYegvKIyHvt4Haf3acuZ/ZK9jiMiIh5QgYegB2dkUFnjuHuc7jQmIhKpVOAhZuGmXbyzNIcJp3Sjc5tmXscRERGPqMBDSHWN4w//XUX7+DiuO62713FERMRDKvAQ8so3W1i9vYg7xvahaazOPxQRiWQq8BCxe28Ff/9gDUO7tWZs/xSv44iIiMdU4CHiHx+uobisinvP7af7fIuIiAo8FKzKKeTl+Vu4bGhneifrZiUiIqICD3rOOf44fTUJTWP5zek9vY4jIiJBQgUe5KYvy+GbTbu49axexDfVjGsiIlJLBR7E9pZX8cCMDI5PjecX6R29jiMiIkFE1yIFsX99kkVeUTlPXjqYqCiduCYiIv+jEXiQ2lBQwqS5G/jpoFQGdWrldRwREQkyKvAg9af3VtO4UTS3je7ldRQREQlCKvAg9HFGHp+uKeCm03vQtkWc13FERCQIqcCDTFllNfe9t5q0ts25/OQuXscREZEgpZPYgsykuRvZvHMfL111EjHR+v1KREQOLuANYWZnm9kaM8sys9sP8vooMys0s6W+rz8EOqNXcvaU8vgnWZzdL5kRPRK9jiMiIkEsoCNwM4sG/g2cAWQDC8xsunNu9QGrznHOjQtktmDwwIwMapzjzrF9vI4iIiJBLtAj8CFAlnNug3OuApgGnBfgDEFp3oadvLd8OxNP7U7H1k29jiMiIkEu0AXeAdha53m2b9mBhpnZMjObaWb9DrYhM5tgZgvNbGFBQUFDZA2Yquoa7p2+ig4JTfjVqO5exxERkRAQ6AI/2HRi7oDni4HOzrkBwL+Adw62Iefc0865dOdcelJSkp9jBtbU+VvIzC3m7nF9iIuJ9jqOiIiEgEAXeDZQd1LvVCCn7grOuSLnXInv8QwgxszC9oyunSXl/OODNYxIS+SsfslexxERkRAR6AJfAPQws65mFgtcBEyvu4KZJZuZ+R4P8WXcGeCcAfP3D9awr6Kae8/ti+9ji4iI/KCAnoXunKsysxuA2UA0MNk5t8rMJvpefwr4GfArM6sCSoGLnHMH7mYPC8uz9zBtwVauHN6VtLYtvI4jIiIhJOATufh2i884YNlTdR4/Djwe6FyBVlPjuHf6Kto0a8yNp/fwOo6IiIQYTfXlkbeXbGPxlj3cdnYvWsbFeB1HRERCjArcA8VllTw4M5MTOiXw00GpXscREZEQpLnQPfDYx+vYubecyb9MJypKJ66JiMiR0wg8wLLyi3nuy01cmN6R41MTvI4jIiIhSgUeQM45/vjuaprERvO7s3p5HUdEREKYCjyAPlidx5x1O/jtGT1JbN7Y6zgiIhLCVOABUlZZzZ/eW02vdi24bGhnr+OIiEiI00lsAfKfzzeQvbuUV64ZSqNo/d4kIiLHRk0SANm79/HEZ1mMPT6FYd3beB1HRETCgAo8AP46aw1mcOeYPl5HERGRMKECb2ArsguZviyHq0d0o31CE6/jiIhImFCBNyDnHA/OzKB1s1iuPbWb13FERCSMqMAb0OdrC/hq/U5+/aM0Wmi+cxER8SMVeAOprnE8NDOTTq2bcslJumxMRET8SwXeQN5eso3M3GJuOasXsY301ywiIv6lZmkAZZXVPPzBGo5PjWds/xSv44iISBhSgTeAF77aRE5hGbeP7q27jYmISINQgfvZnn0V/PvTLE7rlcTJ3RO9jiMiImFKBe5n//40i5LyKm4frUlbRESk4ajA/Wjrrn288NVmfjoolV7JLbyOIyIiYUwF7kcPf7gWM/jtmT29jiIiImFOBe4nK7cV8s7SbVw5oisp8ZoyVUREGpYK3E/+MiuT+CYxTDy1u9dRREQkAqjA/eCLtQXMWbeDG05LI76JpkwVEZGGpwI/RjW+KVNTWzXhsmGaMlVERAJDBX6M/rtsG6u3F3HLWb1o3Cja6zgiIhIhVODHoKyymr/PXstxHVpyzvHtvY4jIiIRRAV+DF78ejPb9pTy+9F9NGWqiIgElAr8KBXuq+TxT7MY2TOJ4WmaMlVERAJLBX6Unvgsi6KySm4/u7fXUUREJAKpwI/Ctj2lPPfVJs4/oQN927f0Oo6IiEQgFfhRePiDtQDcfGYvj5OIiEikUoEfodU5Rby1JJsrTu5ChwRNmSoiIt5QgR+hv8zKpGVcDNeNSvM6ioiIRDAV+BH4MmsHn68t4PrTuhPfVFOmioiId1Tg9VRT43hwZgYdEpowflgXr+OIiEiEU4HX07vLc1i5rYibz+xJXIymTBUREW+pwOuhvKqav81eQ5+UlvxkYAev44iIiKjA6+OleVvI3l3K70f31pSpIiISFFTgP6CwtJLHP1nHiLRERvZM8jqOiIgIoAL/QU99vp7d+yq5fbSmTBURkeChAj+M7YWlTJ67kZ8MbM9xHeK9jiMiIvItFfhhPPzBWpzTlKkiIhJ8VOCHsCa3mDcXZzN+WGc6tm7qdRwREZHvUIEfwl9mZdKscSOuP01TpoqISPAJeIGb2dlmtsbMsszs9oO8bmb2mO/15WY2KNAZv16/k08y87n+tDRaNYsN9NuLiIj8oIAWuJlFA/8GRgN9gYvNrO8Bq40Gevi+JgBPBjLj/ilTU+Lj+OXJXQL51iIiIvUW6BH4ECDLObfBOVcBTAPOO2Cd84AprtY8IMHMUgIV8P0V21meXchvz9CUqSIiErwCXeAdgK11nmf7lh3pOpjZBDNbaGYLCwoK/BZwZ0k5AzsmcMGgVL9tU0RExN8CXeAHm4fUHcU6OOeeds6lO+fSk5L8N0PaL4d35a1fnUy0pkwVEZEgFugCzwY61nmeCuQcxToNSvOdi4hIsAt0gS8AephZVzOLBS4Cph+wznRgvO9s9KFAoXNue4BzioiIBLVGgXwz51yVmd0AzAaigcnOuVVmNtH3+lPADGAMkAXsA64IZEYREZFQENACB3DOzaC2pOsue6rOYwdcH+hcIiIioUQzsYmIiIQgFbiIiEgIUoGLiIiEIBW4iIhICFKBi4iIhCAVuIiISAhSgYuIiIQgFbiIiEgIstp5U0KbmRUAm/24yURghx+3F4zC/TOG++cDfcZwEO6fD/QZ/aGzc+57d+0KiwL3NzNb6JxL9zpHQwr3zxjunw/0GcNBuH8+0GdsSNqFLiIiEoJU4CIiIiFIBX5wT3sdIADC/TOG++cDfcZwEO6fD/QZG4yOgYuIiIQgjcBFRERCkApcREQkBKnAD2BmZ5vZGjPLMrPbvc7jb2Y22czyzWyl11kagpl1NLNPzSzDzFaZ2Y1eZ/I3M4szs2/MbJnvM/7R60wNwcyizWyJmb3ndZaGYGabzGyFmS01s4Ve52kIZpZgZm+YWabv/+QwrzP5i5n18v3b7f8qMrObAppBx8D/x8yigbXAGUA2sAC42Dm32tNgfmRmI4ESYIpz7jiv8/ibmaUAKc65xWbWAlgE/CTM/g0NaOacKzGzGGAucKNzbp7H0fzKzH4LpAMtnXPjvM7jb2a2CUh3zoXtJCdm9gIwxzn3rJnFAk2dc3u8zuVvvu7YBpzknPPnpGKHpRH4dw0BspxzG5xzFcA04DyPM/mVc+4LYJfXORqKc267c26x73ExkAF08DaVf7laJb6nMb6vsPpN3MxSgbHAs15nkaNjZi2BkcAkAOdcRTiWt8+PgfWBLG9QgR+oA7C1zvNswuyHfyQxsy7ACcB8b5P4n2/38lIgH/jQORdun/ER4FagxusgDcgBH5jZIjOb4HWYBtANKACe8x0KedbMmnkdqoFcBLwS6DdVgX+XHWRZWI1sIoWZNQfeBG5yzhV5ncffnHPVzrmBQCowxMzC5nCImY0D8p1zi7zO0sCGO+cGAaOB632Ht8JJI2AQ8KRz7gRgLxCO5xXFAucCrwf6vVXg35UNdKzzPBXI8SiLHCXfceE3ganOube8ztOQfLskPwPO9jiKPw0HzvUdI54G/MjMXvI2kv8553J8f+YDb1N7CC+cZAPZdfYOvUFtoYeb0cBi51xeoN9YBf5dC4AeZtbV91vVRcB0jzPJEfCd4DUJyHDOPex1noZgZklmluB73AQ4Hcj0NpX/OOd+75xLdc51ofb/4CfOuUs9juVXZtbMd5Ilvt3KZwJhdWWIcy4X2GpmvXyLfgyEzcmkdVyMB7vPoXYXh/g456rM7AZgNhANTHbOrfI4ll+Z2SvAKCDRzLKBe5xzk7xN5VfDgcuAFb5jxAB3OOdmeJjJ31KAF3xnvkYBrznnwvJSqzDWDni79vdNGgEvO+dmeRupQfwamOobEG0ArvA4j1+ZWVNqr1q61pP312VkIiIioUe70EVEREKQClxERCQEqcBFRERCkApcREQkBKnARUREQpAKXCSMmdlrZrbLzJIPWB5tZgvNbJ3vWvJDff+9ZuZ8X5v8mOuuOtvN9td2RSKJClwkvN1A7XTATxyw/HfUzop1tXOutB7bGQac78dcz/m2GU7X54sElApcJIz5pun8DXC+mf0cwMx6AvcC/3HOfV7P7cxzzi3xY65tvtufFvhrmyKRRgUuEuacc1OAWcDjZpZE7VSzBcBtR7tNM+vi2/090cweNLNcMys2s5fMrKmZpZnZbDMrMbMsM7vcTx9HRHxU4CKR4VqgKTAPGAFM9NNd2n4PtAcuB/4AXAg8Re3NOd6ndrf7cmpvKdnPD+8nIj6aC10kAjjntpjZ49TezvEtP84Nv945t390PdvMTqF2LvrLnHMvAZjZQmpvt/gzIKzuLSDiJY3ARSKAmbWktlgdcOL+O2H5wcwDnu+/K9rs/Qucc7uBfL57q14ROUYqcJHI8DegFTAWaAs86Kft7j7gecVhlsf56T1FBBW4SNgzs1OBa4C7nHMzgfuB68zsZG+TicixUIGLhDHfJC3PAguAR32L/wKsBJ713adZREKQClwkvN0HdKZ2wpYaAOdcJXAV0Au408NsInIMVOAiYcrM0qmdxOUh59yKuq855/aPyG/X5V0iocmcc15nEJEgZWb3AvcAMYBzzlX7absGRFM7qcyPnXOp/tiuSCTRCFxE6qMSWO/H7d3p2+Z4P25TJKJoBC4ih2Rm7amdaQ2g/MBd8cew3RSgg+9phXNuuT+2KxJJVOAiIiIhSLvQRUREQpAKXEREJASpwEVEREKQClxERCQEqcBFRERC0P8HlLEKpDSh2CcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np                # use arrays\n",
    "import matplotlib.pyplot as plt   # plot\n",
    "%matplotlib inline\n",
    "\n",
    "# Set up the initial conditions, we will assume all magnitudes are in SI units\n",
    "# and angles are in degrees\n",
    "v0 = 10\n",
    "theta = 45\n",
    "\n",
    "# Trigonometric functions take arguments in radians\n",
    "deg2rad = np.pi / 180\n",
    "\n",
    "# The value of gravity is fixed\n",
    "g = 9.81\n",
    "\n",
    "# Set up the times at which we want to evaluate the position of the ball\n",
    "tmin, tmax, nts = 0, 1, 20\n",
    "ts = np.linspace(tmin, tmax, nts)\n",
    "\n",
    "# Compute the positions\n",
    "xs = v0 * np.cos(theta * deg2rad) * ts\n",
    "ys = - g * ts**2 / 2 + v0 * np.sin(theta * deg2rad) * ts\n",
    "zs = 0\n",
    "\n",
    "# Plot them\n",
    "font = 16   # set font size\n",
    "\n",
    "fig = plt.figure(figsize=(7, 5))\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.plot(xs, ys)\n",
    "ax.set_xlabel(r\"X [m]\", fontsize=font), ax.set_ylabel(r\"Y [m]\", fontsize=font)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see how the trajectory is a parabola, in which the ball gains height, then turn over and falls while it advances horizontally at constant speed.\n",
    "\n",
    "In general, it is better to package the code using functions for re-usability. This way we can use the same code to solve many different problems. It makes sense to use a version of the equations prior to performing any substitutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:37:23.600624Z",
     "start_time": "2021-01-31T23:37:23.592939Z"
    }
   },
   "outputs": [],
   "source": [
    "# First, let's code a helper function that gives us the velocity components\n",
    "# given the speed and aim angle\n",
    "def initial_velocity(v0, theta):\n",
    "    '''Returns velocity components given speed and angle with horizontal'''\n",
    "    import numpy as np\n",
    "\n",
    "    deg2rad = np.pi / 180\n",
    "    v0_x = v0 * np.cos(theta * deg2rad)\n",
    "\n",
    "    # START YOUR CODE HERE (1 line of code)>>>>>>>>>>>>>>>>>>\n",
    "    v0_y = v0 * np.sin(theta * deg2rad)\n",
    "    # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> END YOUR CODE HERE\n",
    "\n",
    "    v0_z = 0\n",
    "\n",
    "    return v0_x, v0_y, v0_z\n",
    "\n",
    "# Second, let's create a function that gives the position of the ball at a\n",
    "# given time t\n",
    "\n",
    "\n",
    "def parabolic_positions(t, r0, v0):\n",
    "    '''Returns positions at time t given initial positions r0 and velocities v0'''\n",
    "\n",
    "    # define gravitational constant\n",
    "    g = 9.81\n",
    "\n",
    "    # unpackage inputs\n",
    "    x0 = r0[0]\n",
    "    y0 = r0[1]\n",
    "    z0 = r0[2]\n",
    "    vx0 = v0[0]\n",
    "    vy0 = v0[1]\n",
    "    vz0 = v0[2]\n",
    "\n",
    "    x = vx0 * t + x0\n",
    "\n",
    "    # START YOUR CODE HERE (1 line of code)>>>>>>>>>>>>>>>>>>\n",
    "    y = -g * t**2 / 2 + vy0 * t + y0\n",
    "    # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> END YOUR CODE HERE\n",
    "\n",
    "    z = vz0 * t + z0\n",
    "\n",
    "    return x, y, z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:37:24.712565Z",
     "start_time": "2021-01-31T23:37:24.706126Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7.0710678118654755, 7.071067811865475, 0)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's test the output for the initial velocity helper function\n",
    "initial_velocity(10, 45)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Expected Output**: \n",
    "<table style = \"width:40%\">\n",
    "    <tr>\n",
    "    <td>** initial_velocity(10, 45) **</td> \n",
    "        <td>(7.0710678118654755, 7.071067811865475, 0) </td> \n",
    "    </tr>\n",
    "\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:37:25.541612Z",
     "start_time": "2021-01-31T23:37:25.535316Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.5, 1.27375, 0.0)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's test the output for the computation of parabolic positions\n",
    "parabolic_positions(0.5, (0, 0, 0), (5, 5, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Expected Output**: \n",
    "<table style = \"width:40%\">\n",
    "    <tr>\n",
    "    <td>** parabolic_positions(0.5, (0,0,0), (5,5,0)) **</td> \n",
    "        <td>(2.5, 1.27375, 0.0) </td> \n",
    "    </tr>\n",
    "\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:37:27.478824Z",
     "start_time": "2021-01-31T23:37:27.317936Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAFgCAYAAABEyiulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3iV9f3/8ec7ISHMBEgggbDDFkGICIKIrYuhVjvUn4p1IVX71dY66qjWWrXLqrVqVVBRFLdFZbgVVJA9EyDsEDJYGZCdz++PHGxEwAAn5z7j9biuXJxznzv3eR24yCufe3xuc84hIiIioSXK6wAiIiJy5FTgIiIiIUgFLiIiEoJU4CIiIiFIBS4iIhKCGnkdwB8SExNdly5dvI4hIiLid4sWLdrhnEs6cHlYFHiXLl1YuHCh1zFERET8zsw2H2y5dqGLiIiEIBW4iIhICFKBi4iIhCAVuIiISAhSgYuIiIQgFbiIiEgIUoGLiIiEIBW4iIhICFKBi4iIhCAVuIiISAhSgYuIiISgsJgLXUQkmDjnKK2spryyhorqGiqq6vx5wOPyA1+rqv7u63UeR0cZ7VrG0T4hjpT4JrSPb0K7+MY0bhTt9UcWD6jARUSOknOOvKJy1uYVszavmKz8EtbmFbMuv4Tisiq/vEdsdBSxjWq/KqtrDrrdxOaxpMQ3ITk+jvbxcaQkNCElvrbkU+LjSI6PIyZaO1zDjQpcROQHOOfILSpjXV7JYYu6dbNYerRtznkD25PaqimNfcW7v4QbN4oiJvq7y/Yvj42O/vb5/tdjog0z+06WveVVbC8sY3thae2fe/73ePPOvcxbv5Pi8u+WvBkkNm9cW+7xTUhJiPu24Du1bkrf9i1V8CEooAVuZh2BKUAyUAM87Zx79IB1RgH/BTb6Fr3lnLsvkDlFJDLtL+q1eSWsyytmXV4J6/IPXdQ/GdiBHu2a06NtC3q2a06b5o0bPGOzxo1Ia9uctLbND7lOcVkluYVl5BSWkVtYSk6dks8qKGHOugL2VlT/b5ux0aR3ac3Qbm0Y2q01/TvE00iFHvQCPQKvAm52zi02sxbAIjP70Dm3+oD15jjnxgU4m4hEmILicj7NzGfR5t2szS8mK6/kO6PXNs1iSfMVdc92zUkLYFEfixZxMbSIi6FHuxYHfd05R3F5Fdv3lLEuv5j5G3Yxb8NO/jIrE6gt9BO77i/0NhzXvqUKPQgFtMCdc9uB7b7HxWaWAXQADixwERG/c86RlV/Chxl5fLQ6jyVb9+BcnRH1CaFV1EfLzGgZF0PL5Bh6Jbdg3PHtgdpfaOZv3Mm8DTuZt2EXD82sLfTmjRtxYpdWDO3WhmHd29A3RYUeDMw5580bm3UBvgCOc84V1Vk+CngTyAZygN8551Yd5PsnABMAOnXqNHjz5s0NH1pEQk5VdQ0LN+/mo9V5fJiRx+ad+wDo3yGeM/q24/Q+7eiT0uJ7x5oF8ovLvh2df71hJxsK9gLQonEj3wi9NcO6JdK3fUuio/T311DMbJFzLv17y70ocDNrDnwO/Nk599YBr7UEapxzJWY2BnjUOdfjcNtLT093CxcubLjAIhJSSsqr+HxNAR9l5PFJZj6FpZXERkdxclobTu/Tjh/3aUtKfBOvY4ac/KIy5m3cxdfrdzJ/w0427PhfoQ/p2pph3Wt3ufdJUaH7U9AUuJnFAO8Bs51zD9dj/U1AunNux6HWUYGLSM6eUj7OyOPDjHzmrd9JRXUNCU1j+FHvtpzRpx2n9EyieWNdeONPeUVlvt3ttbvcN/oKvVXTGM4b2IGLhnSkd3JLj1OGvqAocKvdR/UCsMs5d9Mh1kkG8pxzzsyGAG8And1hgqrARSKPc45VOUV8lJHHh6vzWJVTeySuS5umnNG3HWf0TWZQpwQdqw2g3MLaQv8wI48PV+VRUV3DgI4JXHRiR84Z0F6/QB2lYCnwEcAcYAW1l5EB3AF0AnDOPWVmNwC/ovaM9VLgt865rw63XRW4SGSoqKph3oadfOQ7CS2nsAwzGNypFaf7jmd3T2qm49lBYNfeCt5eso1XF2xhbV4JTWOjGXd8Chee2IlBnRL0b3QEgqLAG4oKXCS87Swp56V5W3hx3iZ2lFTQJCaaU3okcnrfdvyod1sSw/Rs8XDgnGPJ1j28+s1W3l2ew76Kanq0bc6FJ3bkgkGptG4W63XEoKcCF5GQsy6vmMlfbuTNxduoqKrhtF5JXHJSZ0b0SCQuRvN/h5qS8ireX57DtAVbWbJlDzHRxpn9krnoxI4M755IlE58OygVuIiEBOccc7N2MGnuRj5bU0DjRlFcMCiVq0Z0Ia3twScmkdCzJreYVxds5a0l2ezZV0lqqyb8Ir0jP09P1RUCB1CBi0hQK6+q5r9Lc5g8dyOZucUkNm/M+GGdueSkTmE7oYrU/rt/sCqPVxdsZW7WDqIMTu2ZxIUnduLHfdpqjnZU4CISpHaWlDN1/hamfL2ZHSXl9E5uwVUjunLuwPa6TWaE2bprH68t3MrrC7PJLSojsXksPx2cyoXpHemWdOi538OdClxEgkpWfjGT5m7ircXZlFfVMKpXEleP6MbwtDY6QznCVdc4vlhbwLQFW/g4I5+qGseQrq258cc9GJ6W6HW8gFOBi4jnnHN8mbWTZ+du0PFtqZf84jLeWryNKV9tIqewjOFpbbj1rN4M6JjgdbSAUYGLiGfKq6qZvjSHSd8e345l/LAuOr4t9VZeVc3UeVt4/NMsdu2tYPRxydx8Zq/D3lY1XKjARSTgdu2tYOq8zbzgO77dq10LrjqlK+cOaK/LwOSolJRXMWnORp6Zs4F9FVX8bHAqN57ekw4J4XvmugpcRAJm994KHv14Ha98s+Xb49tXjejKiLREHd8Wv9hZUs4Tn63nxa9r70R52bDOXDeqe1ju0VGBi0iDq6qu4ZVvtvCPD9dSXFbFzwalcvUpXenRTse3pWFs21PKox+t5Y1F2TSJieaakd24+pRuYTXvugpcRBrU1+t38sd3V5GZW8zJ3dtwzzn96JWs4pbAyMov5h8frGXmylxaN4vl+tPSuOSkTmFxqEYFLiINInv3Ph6ckcn7K7bTIaEJd43tw9nHJWtXuXhi2dY9/G32GuZm7aBDQhNuPL0HF5zQIaTvSqcCFxG/Kqus5qnP1/PkZ+sxg+tGpTFhZLewGPFI6Psyawd/nZXJsuxC0to253dn9uSsfqH5i6UKXET8wjnHzJW5/Pn9DLbtKWXc8Sn8fkyfsD4LWEKTc47Zq3L52+w1rC/Yy4COCdx6Vq+QmwxGBS4ixywzt4h7p69i3oZd9E5uwb3n9mNotzZexxI5rKrqGt5aso1HPlxLTmEZI9ISueWsXiEzGYwKXESO2p59FTz84VpemreZlk1iuPnMXlx8YseQPq4okaesspqp87fw7zqTwdw1rm/Q7z1SgYvIEauucbz8zRb+8cEaikoruXRoZ357Rk8SmsZ6HU3kqBWXVTJp7kae/mIDUWbcNbYPF57YMWiPj6vAReSIzN+wk3vfXU3G9iKGdmvNPef0o09KS69jifjN1l37uOWNZczbsIuRPZN46IL+tA/C0bgKXETqJWdPKQ/MyOC95bWXhd05tg+jdVmYhKmaGsdL8zfz4IxMGkUZd43rwy/Sg2s0rgIXkcMqq6zm6S828MRnWTgHvxrVnWtHdqdJrC4Lk/C3ZWftaHz+xl2c2jOJh37an5T44BiNq8BF5JA+WJXLfe+tJnt3KWP7p/D7Mb1JbdXU61giAVVT43hx3mYemplJo2jj7nF9+fngVM9H44cq8PCZLFZEjti+iirue3c10xZspXdyC165ZijDuuuyMIlMUVHG5Sd3YVSvJG55fTm3vrGcmSu28+AFx5McH+d1vO/RCFwkQq3KKeT/XlnChh17+dWp3fnNGT2J0WVhIkDtaPyFrzfxl1mZxERHcc85/fjpoA6ejMYPNQLX/1aRCOOcY9LcjZz/768oKa9i6lUncevZvVXeInVERRlXDO/KrBtH0ju5Bb97fRlXvbCQvKIyr6N9S/9jRSLIjpJyrnx+AX96bzUjeyYx88aRnBxi00qKBFKXxGa8OmEYd4/ry1frd3DGw5/z1uJsgmHvtQpcJEJ8sbaAsx+Zw5frd/Kn8/rxzPjBtG6mCVlEfkhUlHHViK7MvHEkPdu14LevLeOaKQvJ93g0rgIXCXMVVTU8MCOD8ZO/oXWzGKbfMJzLhnXx/MxakVDTNbEZr147jLvG9mHOuh2c8c8veGfJNs9G4ypwkTC2ccdefvrkVzz9xQYuHdqJ6TeMoHeyZlMTOVrRUcbVp3Rjxo2n0D2pGTe9upQJLy4ivzjwo3EVuEgYcs7x+sKtjH1sDlt37+M/lw3m/p/01726Rfyke1JzXp94MneO6cMXaws4859f8N+lgR2Nq8BFwkxRWSX/N20pt7yxnONT45l54ymc1S/Z61giYSc6yrhmZO1ovGtiM26ctpRrX1xEQXF5QN5fBS4SRhZt3s2YR+cwY8V2bjmrF1OvHho000GKhKvuSc15Y+LJ3DGmN5+tLeCbjbsC8r6aiU0kDFTXOJ74NItHPl5H+4Q4Xp84jEGdWnkdSyRiREcZE0Z255wB7QP2S7MKXCTEbS8s5aZpS5m/cRfnDmjP/ecfR8u4GK9jiUSkQO7xUoGLhLBZK3O57c3lVFbX8PefD/BsqkcRCTwVuEgIKq2o5v73VzN1/hb6d4jnsYtPoGtiM69jiUgAqcBFQkxmbhG/fnkJ6/JLmDCyG787sxexjXQ+qkikUYGLhJBZK3O56dUlNG8cw5QrhzCyZ5LXkUTEIypwkRCw/w5if56RwYDUBJ4Zn05Si8ZexxIRD6nARYJcVXUNf3x3NS/O28zo45L554UDNaOaiKjARYJZSXkVv355MZ+uKeDakd247ezeREXpLHMRUYGLBK3cwjKufH4Ba/KK+fP5x3HJSZ29jiQiQUQFLhKEVucUceXzCyguq2TS5emM6tXW60giEmRU4CJB5tM1+dwwdTEt4mJ4feLJ9G2v23+KyPepwEWCyEvzNnPP9FX0Tm7BpMtPJDk+zutIIhKkVOAiQaCmxvHgzAyembORH/Vuy78uPoFmjfXfU0QOLaDTN5lZRzP71MwyzGyVmd14kHXMzB4zsywzW25mgwKZUSTQSiuquW7qYp6Zs5Hxwzrz9GWDVd4i8oMC/VOiCrjZObfYzFoAi8zsQ+fc6jrrjAZ6+L5OAp70/SkSdgqKy7l6ykKWZ+/h7nF9uXJ4F92MRETqJaAF7pzbDmz3PS42swygA1C3wM8DpjjnHDDPzBLMLMX3vSJhIyu/mF8+t4AdJeU8delgzuqX7HUkEQkhnu2nM7MuwAnA/ANe6gBsrfM827fsOwVuZhOACQCdOnVqqJgiDeKr9TuY+OIiYhtF89q1wzg+NcHrSCISYjy5hZGZNQfeBG5yzhUd+PJBvsV9b4FzTzvn0p1z6UlJuqGDhI43FmUzftI3JMfH8c71J6u8ReSoBHwEbmYx1Jb3VOfcWwdZJRvoWOd5KpATiGwiDck5xz8/XMtjn2QxIi2RJy4dRMu4GK9jiUiICvRZ6AZMAjKccw8fYrXpwHjf2ehDgUId/5ZQV15VzU2vLuWxT7K4ML0jz11xospbRI5JoEfgw4HLgBVmttS37A6gE4Bz7ilgBjAGyAL2AVcEOKOIX+3eW8G1Ly7im027uOWsXlw3qrvONBeRYxbos9DncvBj3HXXccD1gUkk0rA27djLFc8vYNueUv518QmcM6C915FEJExotgiRBpKxvYhLnp2Pc46Xrz6J9C6tvY4kImFEBS7SAPaXd+NGUbx8zVC6JjbzOpKIhBkVuIifZeb+r7xfuWYoXVTeItIAPLkOXCRcZeYW8f+emU9stMpbRBqWClzET+qW97QJKm8RaVgqcBE/UHmLSKCpwEWO0ZrcYpW3iAScClzkGKzJLebiZ+bVHvNWeYtIAKnARY7SgeWtS8VEJJBU4CJHoXa3+Txiok3lLSKeUIGLHKH95d0o2pg2YZjKW0Q8oQIXOQJr81TeIhIcVOAi9bQ2r5iLn64t71c0PaqIeEwFLlIPB5Z3t6TmXkcSkQinAhf5Aft3m0dHqbxFJHiowEUOY395R5kxbYLKW0SChwpc5BDWqbxFJIipwEUOYl1e7SQtKm8RCVYqcJED1C3vV1TeIhKkVOAiddSW9/xvy7u7yltEgpQKXMRnf3mbofIWkaCnAhcBNu/c+215T1N5i0gIUIFLxNu1t4JfPreAqpoaXrlG5S0ioaGR1wFEvFRWWc01UxaybU8pL199EmltVd4iEho0ApeIVVPj+O1rS1m8ZTePXjiQ9C6tvY4kIlJvKnCJWA/MyGDGilzuHNOH0f1TvI4jInJEVOASkZ77ciPPzt3IFcO7cPUp3byOIyJyxFTgEnFmrczlvvdWc1a/dtw1tq/XcUREjooKXCLKos27uXHaEgZ2TODRi04gOsq8jiQiclRU4BIxNu7Yy9UvLCAlPo5nx6cTFxPtdSQRkaOmApeIsLOknCue+wYz4/krhtCmeWOvI4mIHBMVuIS90opqrp6ykO2FZTx7eTpdEpt5HUlE5JhpIhcJa9U1jpteXcLSrXt48pLBDOrUyutIIiJ+oRG4hC3nHH96bzWzV+Xxh3F9Ofu4ZK8jiYj4jQpcwtakuRt5/qtNXD2iK1cM7+p1HBERv1KBS1iasWI7f56RwZj+ydwxpo/XcURE/E4FLmFn4aZd3PTqUgZ3asXDvxhIlK71FpEwpAKXsLKhoISrpywkNaEJz+habxEJYypwCRs7Ssr55XMLiPZd692qWazXkUREGowuI5OwsK+iiqueX0B+cRnTJgyjU5umXkcSEWlQGoFLyKuucfzfK0tZsa2Qf108iIEdE7yOJCLS4DQCl5DmnOPe6av4KCOP+87rxxl923kdSUQkIDQCl5D2zJwNvDhvM9eO7Mb4YV28jiMiEjAqcAlZ7y7L4YEZmYw7PoXbzu7tdRwRkYBSgUtI+mbjLm5+bRlDurTm7z8foGu9RSTiBLTAzWyymeWb2cpDvD7KzArNbKnv6w+BzCehYdOOvVwzZSEdWzfh6fGDda23iESkQJ/E9jzwODDlMOvMcc6NC0wcCTWlFdVMfGkRZvD8FUNIaKprvUUkMgV0BO6c+wLYFcj3lPDhnOOOt1ewJq+Yxy46gY6tda23iESuYDwGPszMlpnZTDPrd6iVzGyCmS00s4UFBQWBzCceeXHeZt5eso2bz+jJyJ5JXscREfFUsBX4YqCzc24A8C/gnUOt6Jx72jmX7pxLT0rSD/Nwt2jzLu57dzWn92nLdaPSvI4jIuK5oCpw51yRc67E93gGEGNmiR7HEo/lF5dx3dTFdGjVhH/o7mIiIkCQFbiZJZuZ+R4PoTbfTm9TiZcqq2u44eUlFJZW8tSlg4lvEuN1JBGRoBDQs9DN7BVgFJBoZtnAPUAMgHPuKeBnwK/MrAooBS5yzrlAZpTg8tdZmXyzcRf/vHAAfVJaeh1HRCRoBLTAnXMX/8Drj1N7mZkI7y/fzjNzNnL5sM6cf0Kq13FERILKIQvczA53rfbh3O2c23yU3ysCQFZ+Mbe8sYxBnRK4c2xfr+OIiASdw43ALwVygfIj2F4n4BFABS5HrbiskgkvLqJpbDRPXDKY2EZBdaqGiEhQ+KFd6D9xzn1Tnw2ZWSOg4tgjSSRzznHrG8vZvHMfL111EsnxcV5HEhEJSocb2nwMFB7Btqp931N0TIkkoj0zZwMzV+Zy+9m9Gda9jddxRESC1iFH4M65M45kQ76zxY/oe0Tq+mr9Dh6amcmY/slcfUpXr+OIiAQ1HVyUoLC9sJRfv7yEronN+OvPBuCbDkBERA6h3peR+Y5xDwM6At87MOmcm+zHXBJBKqpquG7qYsoqq/nPZUNp3jjQN8kTEQk99fpJaWaDgLeBVOBgQyMHqMDlqNz//mqWbNnDE5cMIq1tC6/jiIiEhPoOdZ4CSoCfAJnobHPxk7cWZzPl681MGNmNMf1TvI4jIhIy6lvgfYFf+G4wIuIXq3OKuOPtFQzt1ppbz+rldRwRkZBS35PY1gLNGjKIRJbCfZVMfGkR8U1i+NfFg2gUrfMpRUSORH1/at4B3GVmnRoyjESGmhrHb15byvbCUp64ZDBJLRp7HUlEJOTUaxe6c26WmY0C1pnZWmD391dxp/o7nISnxz/N4pPMfO47rx+DO7fyOo6ISEiq71notwO3AgXUzrRW3ZChJHx9tiaff360lvNP6MBlQzt7HUdEJGTV9yS2m4D/ADc451TeclS27trHjdOW0qtdCx44v78maxEROQb1PQbeFHhd5S1Hq6yyml9NXUSNczx16WCaxEZ7HUlEJKTVt8BnUjsLm8gRc85x9zsrWbmtiH/+YiBdEnVBg4jIsarvLvRHgOd9uzxn8f2T2HDObfBjLgkj0xZs5fVF2fz6R2mc3red13FERMJCfQv8S9+ffwLuO8Q62icq37NyWyH3/HcVp/RI5KbTe3odR0QkbNS3wK+kdr5zkXorrajmxmlLaNUshkcvOoHoKJ20JiLiL/W9Dvz5Bs4hYejBmRmsL9jLi1cNoXWzWK/jiIiEFc1fKQ3i08x8pny9matGdOWUHklexxERCTuHLHAzm2xmXeu7Ias1WdOtyo6Scm55Yzm9k1twi25SIiLSIA43Ar8cSDzCbR3p90iYcc5x+5srKCqt5JGLBhIXo3MbRUQawuGOgRvwpJkV1XNbOkNJmLZgKx9l5HHX2D70Tm7pdRwRkbB1uAL/gtozz4+kmL8Aio8pkYSsDQUl3PfuakakJXLl8HoffRERkaNwyAJ3zo0KYA4JcZXVNfzm1aXENori7z8fQJQuGRMRaVD1vQ5c5LAe+3gdy7ILeeKSQSTHx3kdR0Qk7OkyMjlmCzft4t+fZvGzwamM6Z/idRwRkYigApdjUlxWyW9eW0qHVk2455y+XscREYkY2oUux+Te6avZtruU1ycOo0VcjNdxREQihkbgctTeX76dNxdnc8NpaQzu3NrrOCIiEeVwM7H9OpBBJLTkFpZxx9srGNAxgV//uIfXcUREIs7hRuCPmNkXZpYWsDQSEmpqHDe/vpSKqhoeuXAgMdHakSMiEmiH+8n7YyAFWGZmvzMzXdgrAEz+ciNfZu3kD+f0pWtiM6/jiIhEpEMWuHPuM6A/8DjwAPC1mek04wiXsb2Iv85awxl923HRiR29jiMiErEOexa6c64MuM3MXgWeBRab2TtAxfdXdZc3UEYJEmWV1dw0bSktm8Tw0AX90U4ZERHv1PcysixgKTAQOIWDFLg/Q0lw+tvsNazJK+a5K06kTfPGXscREYloP1jgZnYu8ATQFLjGOTepwVNJ0JmzroBJczdy+bDOnNarrddxREQi3uEuI0vy7Tp/G1gMHKfyjky791bwu9eXkda2Ob8f08frOCIiwuFH4JlADTDeOTc1QHkkyDjnuOPtFezaW8Gky08kLiba60giIsLhLyP7BOin8o5sbyzKZubKXG4+sxfHdYj3Oo6IiPgc7n7gPw9kEAk+W3bu497pqxjarTXXnNLN6zgiIlKHptCSg6qqruE3ry0lKsr4xy8GEh2lS8ZERIKJ7kYmB/XEZ+tZtHk3j140kA4JTbyOIyIiBwjoCNzMJptZvpmtPMTrZmaPmVmWmS03s0GBzCe1lm7dw6Mfr+MnA9tz3sAOXscREZGDCPQu9OeBsw/z+migh+9rAvBkADJJHXvLq7hp2hKSW8bxx/OO8zqOiIgcQkAL3Dn3BbDrMKucB0xxteYBCWaWEph0AnD/+6vZvGsfD/9iAPFNYryOIyIihxBsJ7F1ALbWeZ7tW/Y9ZjbBzBaa2cKCgoKAhAt3n2Tm8co3W5l4andO6tbG6zgiInIYwVbgBzvV+aDzrDvnnnbOpTvn0pOSkho4VvgrLqvkzrdX0qtdC35zek+v44iIyA8ItrPQs4G696hMBXI8yhJR/jZ7DblFZTx56WBiGwXb73UiInKgYPtJPR0Y7zsbfShQ6Jzb7nWocLdw0y5enLeZK07uysCOCV7HERGRegjoCNzMXgFGAYlmlg3cA8QAOOeeAmYAY6i9fek+4IpA5otE5VXV3PbmctrHN+HmM7XrXEQkVAS0wJ1zF//A6w64PkBxBPj3p+tZX7CXF64cQrPGwXZERUREDiXYdqFLAK3JLebJz7K44IQOnNpTJwKKiIQSFXiEqq5x3PbmclrExXDXuL5exxERkSOkAo9QU77exNKte7jnnL60bhbrdRwRETlCKvAIlL17H3+bvYbTeiVx7oD2XscREZGjoAKPMM457ny79l4y95/fHzPdJlREJBSpwCPMf5fm8PnaAm49q5duEyoiEsJU4BFkZ0k5f3x3FSd0SuCyYV28jiMiIsdABR5B/vTeakrKq/jLT48nOkq7zkVEQpkKPEJ8uiafd5bmcN2oNHq2a+F1HBEROUYq8Aiwt7yKu95eSVrb5lx3Wnev44iIiB9o7swI8LfZa8gpLOWNiSfTuFG013FERMQPNAIPc4u37OaFrzcxfmhnBndu5XUcERHxExV4GKuoquH2N5eT3DKOW87u7XUcERHxI+1CD2NPfb6etXklTP5lOs11pzERkbCiEXiYysov5vFPsjh3QHt+1Lud13FERMTPVOBhqKbGcdubK2jaOJo/nKM7jYmIhCMVeBiaOn8zizbv5u6xfUls3tjrOCIi0gBU4GEmZ08pD83M5JQeiVwwqIPXcUREpIGowMOIc46731lJjYMHdKcxEZGwpgIPI+8t387HmfncfGZPOrZu6nUcERFpQCrwMLF7bwX3Tl/FgNR4rhje1es4IiLSwHRxcJi4//0MCksreenqk3SnMRGRCKAReBiYs66ANxdnM/HU7vRJael1HBERCQAVeIjbV1HFHW+voFtSM274UZrXcUREJEC0Cz3EPfzBWrbuKuW1a4cRF6M7jYmIRAqNwEPYsq17mPzlRi45qRNDurb2Oo6IiASQCjxEVVbXcNuby0lq0ZjbRutOYyIikUa70EPU5Lkbycwt5unLBtMyLsbrOCIiEiUZXg0AAA8VSURBVGAagYegvKIyHvt4Haf3acuZ/ZK9jiMiIh5QgYegB2dkUFnjuHuc7jQmIhKpVOAhZuGmXbyzNIcJp3Sjc5tmXscRERGPqMBDSHWN4w//XUX7+DiuO62713FERMRDKvAQ8so3W1i9vYg7xvahaazOPxQRiWQq8BCxe28Ff/9gDUO7tWZs/xSv44iIiMdU4CHiHx+uobisinvP7af7fIuIiAo8FKzKKeTl+Vu4bGhneifrZiUiIqICD3rOOf44fTUJTWP5zek9vY4jIiJBQgUe5KYvy+GbTbu49axexDfVjGsiIlJLBR7E9pZX8cCMDI5PjecX6R29jiMiIkFE1yIFsX99kkVeUTlPXjqYqCiduCYiIv+jEXiQ2lBQwqS5G/jpoFQGdWrldRwREQkyKvAg9af3VtO4UTS3je7ldRQREQlCKvAg9HFGHp+uKeCm03vQtkWc13FERCQIqcCDTFllNfe9t5q0ts25/OQuXscREZEgpZPYgsykuRvZvHMfL111EjHR+v1KREQOLuANYWZnm9kaM8sys9sP8vooMys0s6W+rz8EOqNXcvaU8vgnWZzdL5kRPRK9jiMiIkEsoCNwM4sG/g2cAWQDC8xsunNu9QGrznHOjQtktmDwwIwMapzjzrF9vI4iIiJBLtAj8CFAlnNug3OuApgGnBfgDEFp3oadvLd8OxNP7U7H1k29jiMiIkEu0AXeAdha53m2b9mBhpnZMjObaWb9DrYhM5tgZgvNbGFBQUFDZA2Yquoa7p2+ig4JTfjVqO5exxERkRAQ6AI/2HRi7oDni4HOzrkBwL+Adw62Iefc0865dOdcelJSkp9jBtbU+VvIzC3m7nF9iIuJ9jqOiIiEgEAXeDZQd1LvVCCn7grOuSLnXInv8QwgxszC9oyunSXl/OODNYxIS+SsfslexxERkRAR6AJfAPQws65mFgtcBEyvu4KZJZuZ+R4P8WXcGeCcAfP3D9awr6Kae8/ti+9ji4iI/KCAnoXunKsysxuA2UA0MNk5t8rMJvpefwr4GfArM6sCSoGLnHMH7mYPC8uz9zBtwVauHN6VtLYtvI4jIiIhJOATufh2i884YNlTdR4/Djwe6FyBVlPjuHf6Kto0a8yNp/fwOo6IiIQYTfXlkbeXbGPxlj3cdnYvWsbFeB1HRERCjArcA8VllTw4M5MTOiXw00GpXscREZEQpLnQPfDYx+vYubecyb9MJypKJ66JiMiR0wg8wLLyi3nuy01cmN6R41MTvI4jIiIhSgUeQM45/vjuaprERvO7s3p5HUdEREKYCjyAPlidx5x1O/jtGT1JbN7Y6zgiIhLCVOABUlZZzZ/eW02vdi24bGhnr+OIiEiI00lsAfKfzzeQvbuUV64ZSqNo/d4kIiLHRk0SANm79/HEZ1mMPT6FYd3beB1HRETCgAo8AP46aw1mcOeYPl5HERGRMKECb2ArsguZviyHq0d0o31CE6/jiIhImFCBNyDnHA/OzKB1s1iuPbWb13FERCSMqMAb0OdrC/hq/U5+/aM0Wmi+cxER8SMVeAOprnE8NDOTTq2bcslJumxMRET8SwXeQN5eso3M3GJuOasXsY301ywiIv6lZmkAZZXVPPzBGo5PjWds/xSv44iISBhSgTeAF77aRE5hGbeP7q27jYmISINQgfvZnn0V/PvTLE7rlcTJ3RO9jiMiImFKBe5n//40i5LyKm4frUlbRESk4ajA/Wjrrn288NVmfjoolV7JLbyOIyIiYUwF7kcPf7gWM/jtmT29jiIiImFOBe4nK7cV8s7SbVw5oisp8ZoyVUREGpYK3E/+MiuT+CYxTDy1u9dRREQkAqjA/eCLtQXMWbeDG05LI76JpkwVEZGGpwI/RjW+KVNTWzXhsmGaMlVERAJDBX6M/rtsG6u3F3HLWb1o3Cja6zgiIhIhVODHoKyymr/PXstxHVpyzvHtvY4jIiIRRAV+DF78ejPb9pTy+9F9NGWqiIgElAr8KBXuq+TxT7MY2TOJ4WmaMlVERAJLBX6Unvgsi6KySm4/u7fXUUREJAKpwI/Ctj2lPPfVJs4/oQN927f0Oo6IiEQgFfhRePiDtQDcfGYvj5OIiEikUoEfodU5Rby1JJsrTu5ChwRNmSoiIt5QgR+hv8zKpGVcDNeNSvM6ioiIRDAV+BH4MmsHn68t4PrTuhPfVFOmioiId1Tg9VRT43hwZgYdEpowflgXr+OIiEiEU4HX07vLc1i5rYibz+xJXIymTBUREW+pwOuhvKqav81eQ5+UlvxkYAev44iIiKjA6+OleVvI3l3K70f31pSpIiISFFTgP6CwtJLHP1nHiLRERvZM8jqOiIgIoAL/QU99vp7d+yq5fbSmTBURkeChAj+M7YWlTJ67kZ8MbM9xHeK9jiMiIvItFfhhPPzBWpzTlKkiIhJ8VOCHsCa3mDcXZzN+WGc6tm7qdRwREZHvUIEfwl9mZdKscSOuP01TpoqISPAJeIGb2dlmtsbMsszs9oO8bmb2mO/15WY2KNAZv16/k08y87n+tDRaNYsN9NuLiIj8oIAWuJlFA/8GRgN9gYvNrO8Bq40Gevi+JgBPBjLj/ilTU+Lj+OXJXQL51iIiIvUW6BH4ECDLObfBOVcBTAPOO2Cd84AprtY8IMHMUgIV8P0V21meXchvz9CUqSIiErwCXeAdgK11nmf7lh3pOpjZBDNbaGYLCwoK/BZwZ0k5AzsmcMGgVL9tU0RExN8CXeAHm4fUHcU6OOeeds6lO+fSk5L8N0PaL4d35a1fnUy0pkwVEZEgFugCzwY61nmeCuQcxToNSvOdi4hIsAt0gS8AephZVzOLBS4Cph+wznRgvO9s9KFAoXNue4BzioiIBLVGgXwz51yVmd0AzAaigcnOuVVmNtH3+lPADGAMkAXsA64IZEYREZFQENACB3DOzaC2pOsue6rOYwdcH+hcIiIioUQzsYmIiIQgFbiIiEgIUoGLiIiEIBW4iIhICFKBi4iIhCAVuIiISAhSgYuIiIQgFbiIiEgIstp5U0KbmRUAm/24yURghx+3F4zC/TOG++cDfcZwEO6fD/QZ/aGzc+57d+0KiwL3NzNb6JxL9zpHQwr3zxjunw/0GcNBuH8+0GdsSNqFLiIiEoJU4CIiIiFIBX5wT3sdIADC/TOG++cDfcZwEO6fD/QZG4yOgYuIiIQgjcBFRERCkApcREQkBKnAD2BmZ5vZGjPLMrPbvc7jb2Y22czyzWyl11kagpl1NLNPzSzDzFaZ2Y1eZ/I3M4szs2/MbJnvM/7R60wNwcyizWyJmb3ndZaGYGabzGyFmS01s4Ve52kIZpZgZm+YWabv/+QwrzP5i5n18v3b7f8qMrObAppBx8D/x8yigbXAGUA2sAC42Dm32tNgfmRmI4ESYIpz7jiv8/ibmaUAKc65xWbWAlgE/CTM/g0NaOacKzGzGGAucKNzbp7H0fzKzH4LpAMtnXPjvM7jb2a2CUh3zoXtJCdm9gIwxzn3rJnFAk2dc3u8zuVvvu7YBpzknPPnpGKHpRH4dw0BspxzG5xzFcA04DyPM/mVc+4LYJfXORqKc267c26x73ExkAF08DaVf7laJb6nMb6vsPpN3MxSgbHAs15nkaNjZi2BkcAkAOdcRTiWt8+PgfWBLG9QgR+oA7C1zvNswuyHfyQxsy7ACcB8b5P4n2/38lIgH/jQORdun/ER4FagxusgDcgBH5jZIjOb4HWYBtANKACe8x0KedbMmnkdqoFcBLwS6DdVgX+XHWRZWI1sIoWZNQfeBG5yzhV5ncffnHPVzrmBQCowxMzC5nCImY0D8p1zi7zO0sCGO+cGAaOB632Ht8JJI2AQ8KRz7gRgLxCO5xXFAucCrwf6vVXg35UNdKzzPBXI8SiLHCXfceE3ganOube8ztOQfLskPwPO9jiKPw0HzvUdI54G/MjMXvI2kv8553J8f+YDb1N7CC+cZAPZdfYOvUFtoYeb0cBi51xeoN9YBf5dC4AeZtbV91vVRcB0jzPJEfCd4DUJyHDOPex1noZgZklmluB73AQ4Hcj0NpX/OOd+75xLdc51ofb/4CfOuUs9juVXZtbMd5Ilvt3KZwJhdWWIcy4X2GpmvXyLfgyEzcmkdVyMB7vPoXYXh/g456rM7AZgNhANTHbOrfI4ll+Z2SvAKCDRzLKBe5xzk7xN5VfDgcuAFb5jxAB3OOdmeJjJ31KAF3xnvkYBrznnwvJSqzDWDni79vdNGgEvO+dmeRupQfwamOobEG0ArvA4j1+ZWVNqr1q61pP312VkIiIioUe70EVEREKQClxERCQEqcBFRERCkApcREQkBKnARUREQpAKXCSMmdlrZrbLzJIPWB5tZgvNbJ3vWvJDff+9ZuZ8X5v8mOuuOtvN9td2RSKJClwkvN1A7XTATxyw/HfUzop1tXOutB7bGQac78dcz/m2GU7X54sElApcJIz5pun8DXC+mf0cwMx6AvcC/3HOfV7P7cxzzi3xY65tvtufFvhrmyKRRgUuEuacc1OAWcDjZpZE7VSzBcBtR7tNM+vi2/090cweNLNcMys2s5fMrKmZpZnZbDMrMbMsM7vcTx9HRHxU4CKR4VqgKTAPGAFM9NNd2n4PtAcuB/4AXAg8Re3NOd6ndrf7cmpvKdnPD+8nIj6aC10kAjjntpjZ49TezvEtP84Nv945t390PdvMTqF2LvrLnHMvAZjZQmpvt/gzIKzuLSDiJY3ARSKAmbWktlgdcOL+O2H5wcwDnu+/K9rs/Qucc7uBfL57q14ROUYqcJHI8DegFTAWaAs86Kft7j7gecVhlsf56T1FBBW4SNgzs1OBa4C7nHMzgfuB68zsZG+TicixUIGLhDHfJC3PAguAR32L/wKsBJ713adZREKQClwkvN0HdKZ2wpYaAOdcJXAV0Au408NsInIMVOAiYcrM0qmdxOUh59yKuq855/aPyG/X5V0iocmcc15nEJEgZWb3AvcAMYBzzlX7absGRFM7qcyPnXOp/tiuSCTRCFxE6qMSWO/H7d3p2+Z4P25TJKJoBC4ih2Rm7amdaQ2g/MBd8cew3RSgg+9phXNuuT+2KxJJVOAiIiIhSLvQRUREQpAKXEREJASpwEVEREKQClxERCQEqcBFRERC0P8HlLEKpDSh2CcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's verify that we get the same result we got before with the\n",
    "# direct calculation\n",
    "\n",
    "# Now we can initialize the positions and velocities as\n",
    "r_t0 = (0, 0, 0)\n",
    "v_t0 = initial_velocity(10, 45)\n",
    "\n",
    "# Determine the same times as before\n",
    "\n",
    "# And get the coordinates with a call to the function\n",
    "xs, ys, zs = parabolic_positions(ts, r_t0, v_t0)\n",
    "\n",
    "# Plot result\n",
    "font = 16   # set font size\n",
    "\n",
    "fig = plt.figure(figsize=(7, 5))\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.plot(xs, ys)\n",
    "ax.set_xlabel(r\"X [m]\", fontsize=font), ax.set_ylabel(r\"Y [m]\", fontsize=font)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The question is which launch angle we need to maximize the height. We can see that the solution in this case can be found analytically. There are two times at which the ball is on the floor, the initial one, $t_0=0$ and the final one:\n",
    "\n",
    "$$y(t_f) = 0 \\rightarrow t_f = \\frac{2 v_0 \\sin \\theta}{g}$$\n",
    "\n",
    "The higher the initial vertical velocity, the longer the ball stays in the air, and the stronger the gravity, the less time the ball stays in the air.\n",
    "\n",
    "The maximum reach is then the $x$ position at that time:\n",
    "\n",
    "$$x_{max} = \\frac{2 v_0^2 \\sin \\theta \\cos \\theta}{g} = \\frac{v_0^2 \\sin(2\\theta)}{g}$$\n",
    "\n",
    "The maximum of this function is achieved when:\n",
    "\n",
    "$$\\sin(2\\theta) = 1 \\rightarrow \\theta = \\frac{\\pi}{4} \\equiv 45 \\rm deg$$\n",
    "\n",
    "And the value for that reach is:\n",
    "\n",
    "$$x_{max} = \\frac{v_0^2}{g}$$\n",
    "\n",
    "Check that the units make sense (it has units of lenght)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:37:28.256661Z",
     "start_time": "2021-01-31T23:37:28.251600Z"
    }
   },
   "outputs": [],
   "source": [
    "# We can validate it using the functions we had before\n",
    "\n",
    "# Let's add first another helper function that gives the time-on-air\n",
    "\n",
    "def flight_time(v_0, theta):\n",
    "    import numpy as np\n",
    "\n",
    "    deg2rad = np.pi / 180\n",
    "    g = 9.81\n",
    "\n",
    "    # START YOUR CODE HERE (1 line of code)>>>>>>>>>>>>>>>>>>\n",
    "    t_f = 2 * v_0 * np.sin(theta * deg2rad) / g\n",
    "    # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> END YOUR CODE HERE\n",
    "\n",
    "    return t_f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:37:29.484696Z",
     "start_time": "2021-01-31T23:37:29.480472Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4416040391163047"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's test the output for the flight_time helper function\n",
    "flight_time(10, 45)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Expected Output**: \n",
    "<table style = \"width:40%\">\n",
    "    <tr>\n",
    "    <td>** flight_time(10, 45) **</td> \n",
    "        <td>1.4416040391163047 </td> \n",
    "    </tr>\n",
    "\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:37:30.966891Z",
     "start_time": "2021-01-31T23:37:30.797583Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAFgCAYAAABEyiulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3xT9f7H8ddpmnTvvcumzAooQ1GGIiIIbkUFcSCoV9Crwr1cB070/sS9BRcqKoiAAxcggotVZMiG7r2brjT5/v6o5aqsliY5Sft5Ph55qMkZb46ln5zv+Q5NKYUQQggh3IuH3gGEEEII0XJSwIUQQgg3JAVcCCGEcENSwIUQQgg3JAVcCCGEcEOeegc4mfDwcJWcnGy345nNZvz8/Ox2vPZGrl/ryTVsHbl+rSPXr3X0uH6bN28uUkpF/P19ly/gycnJbNq0yW7HW7t2LcOGDbPb8dobuX6tJ9ewdeT6tY5cv9bR4/ppmpZ+rPelCV0IIYRwQ1LAhRBCCDckBVwIIYRwQy7/DFwIIUTbYrFYyMrKora2Vu8oLRYUFMTvv//ukGN7e3sTHx+P0Whs1vZSwIUQQjhVVlYWAQEBJCcno2ma3nFapLKykoCAALsfVylFcXExWVlZdOjQoVn7SBO6EEIIp6qtrSUsLMztircjaZpGWFhYi1olnH4HrmnaYaASsAINSqkBzs4ghBBCX1K8j9bSa6JXE/pwpVSRTucWQggh3J40oQshhGh3brjhBiIjI+nVq9df3i8pKeG8886jS5cunHfeeZSWlp70WNdffz1LlixxVNTj0pRSzj2hph0CSgEFvKqUeu0Y20wFpgJERUX1X7x4sd3OX1VVhb+/v92O197I9Ws9uYatI9evdVzh+gUFBdG5c2ddM2zYsAE/Pz9uueUWfvnllyPv33fffYSEhHDXXXcxf/58ysrKeOihh458brVaMRgMfznWtGnTGD16NBMmTGh1rv3791NeXv6X94YPH775mI+blVJOfQGxf/wzEtgGnH2i7fv376/sac2aNXY9Xnsj16/15Bq2jly/1nGF67dr1y69IyillDp06JDq2bPnX97r2rWrysnJUUoplZOTo7p27fqXzysqKpTNZlO33XabSklJUWPGjFEXXHCB+vjjj5VSSm3atEmdffbZql+/fmrUqFFHjvXrr7+q3r17q0GDBqm77777qPM2Oda1ATapY9RHpz8DV0rl/PHPAk3TlgFnAOucnaO9a7BYydpdSmFGJVWldQD4h3gRkRhAfPcQPI2GkxxBCCFab+7KnezKqbDrMXvEBvLAuJ6ntG9+fj4xMTEAxMTEUFBQcNQ2y5YtY8+ePWzfvp38/Hx69OjBDTfcgMVi4R//+AfLly8nIiKCDz/8kDlz5rBw4UKmTJnCa6+9xpAhQ5g9e3ar/nxNnFrANU3zAzyUUpV//Pso4KGT7CbsyFxex9avMti1IQdLnRU08AkwAVBTUQ+A0dtAj7NiOe28RPyCvPSMK4QQLmfdunVcffXVGAwGYmNjGTFiBAB79uxhx44dnHfeeUBjc3tMTAxlZWVUVlYyZMgQACZOnMhnn33W6hzOvgOPApb90VXeE3hfKbXKyRnaJaUUO77P5qdPD9BQb6PL6ZF0GxhNTOdgjKbGu21LvZXcfWXs/jmP31Zn8fv6HAZf0pmeQ2NlyIcQwiFO9U7ZUaKiosjNzSUmJobc3FwiIyOPud2xficqpejZsyc//fTTX95vTke4U+HUXuhKqYNKqb5/vHoqpR515vnbq/raBr56fQfrFu8lukMgEx8cyHlTepLYI+xI8QYwmgwk9gxj1I09mfjAQCKTA/n+/T18/cbOxrt1IYRo4y666CLefvttAN5++23Gjx9/1DZnn302ixcvxmq1kpuby5o1awDo1q0bhYWFRwq4xWJh586dhISEEBAQwM8//wyAvTpmyzCyNq7WbGH501s5uLWQIZd0ZtwdqQRH+p50v+AoXy6akcrgiztxYEsBy5/ZSq3Z4oTEQgjheFdffTWDBw9mz549xMfHs2DBAgBmz57NN998Q5cuXfjmm2+O+bz64osvpkuXLvTu3Zvp06dzzjnnAGAymViyZAmzZs2ib9++pKam8uOPPwKwYMECpk6dyuDBg1FKERQU1Oo/g8yF3obVVVv49OmtlOVVc8H0PnToE96i/TVNo9/5SQRH+fLVGztY/sxWwk937rBDIYRwhA8++OCY74eFhfHdd9+dcF9N03jhhReO+Vlqairr1h3dL7tnz5789ttvAMybN48BA1o/CancgbdRDRYrX7y8ndJcM2Nu7d3i4v1nHVMjuOCW3pRkm8lcr7BabHZMKoQQbd/nn39OamoqvXr14ocffuA///lPq48pBbwNUkrx/ft7yNlXxsjrU0jsEdbqYyb3DmfE5BTMBbBu8R47pBRCiPbjyiuvJC0tjR07dvD5558TERHR6mNKAW+Ddv6Qw+6f8hhwYTJdT4+223G7DYwmvAfs2pDLrvU5djuuEEKIlpMC3sYU51Sx/uN9JPYM5YwLm7embEtE9tJISAnhhw/3UppntvvxhRBCNI8U8DbE2mDj2zd3YfI2MHJyDzQP+4/d1jw0Rl7fA0+TgW8W7sJqlefhQgihByngbcjWr9Mpyqxi2DXd8Q00Oew8fkFeDLumG4UZlaR9k+Gw8wghhDg+KeBtRFl+NZu+SKdTv0g6pra+c8TJNJ1n4+eHKS+sdvj5hBDCXmpraznjjDPo27cvPXv25IEHHjjymTstJyoFvA1QSrHuw70YjB4MvbKL08479MqueBg0fvhon9POKYQQreXl5cXq1avZtm0baWlprFq16sgsafPmzWPkyJHs27ePkSNHMm/ePJ3THp8U8DYgfUcxmbtKOGNsB6cuPuIf4sXpYzqQvr2YjJ3FTjuvEEK0hqZpR9ZEt1gsWCyWI3ObL1++nMmTJwMwefJkPv3006P2V0px++2306NHDy688MK/rFi2efNmzjnnHPr378/5559Pbm4uABs3bqRPnz4MHjyYe+65h169erX6zyEzsbk5q9XGhiX7CY7ypdc5cU4/f5/h8ez4IZv1H+/jqu4heBjkO6EQogW+nA152+17zOjecMGJ75ytViv9+/dn//793HbbbQwcOBBwr+VE5betm9v9Yy5l+dUMuaQTBk/n/+80GD0YckknSvOq2fNLntPPL4QQp8JgMJCWlkZWVha//vorO3bsaPa+zVlONDU1lUceeYSsrKxjLidqD3IH7sYaLFY2fn6Y6I6BJLdiqtTW6pgaQWRSAL9+doiup0djMMr3QiFEM53kTtnRgoODGTZsGKtWraJXr16ynKhwjp3rcjCX1TFwfCdd1+vWNI1BEzpRVVLHTpmhTQjh4goLCykrKwOgpqaGb7/9lu7duwPutZyo3IG7qQaLla1fpxPXNZj4biF6xyG+ewgxnYPY+nU6PYfG6tKcL4QQzZGbm8vkyZOxWq3YbDauuOIKxo4dCzQuJ3rFFVewYMECEhMT+fjjj4/a/+KLL2b16tX07t2brl27HrWc6B133EF5eTkNDQ3MnDmTnj17smDBAm6++Wb8/PwYNmyYLCfanu3+KQ9zeT0jp/TQOwrQeBc+4IJkVj6/jT0/59HjrFi9IwkhxDH16dOHrVu3HvMzWU5UOJTNamPr1+lEdQh0ibvvJgk9QolMCmDzV+nYbLJuuBBCNJHlRAUAB9OKqCiqpd+oJF2fff+dpmmcNiqJisIaDm8r0juOEEK4DFlOVKCUIu3bDAIjfEjuq1/P8+PpmBpOQJg3ad/KHOlCCOFIUsDdTN7BCvIPVZA6MgEPB6w21loeBg/6jkwg90A5+Ycq9I4jhBBtlhRwN7N9bRYmH0+6D47RO8pxpQyJwehtYPvaLL2jCCFEmyUF3I2Yy+s4sKWA7oOjMXoZ9I5zXCZvT7oPjGbf5nxqKuv1jiOEEG2SFHA3smt9Djarovc58XpHOalew+KxNSh2bZCJXYQQrslqtXLaaacdGQMO8OCDDxIXF0dqaiqpqal88cUXJz2OLCcqTshmU+xan0N89xCCo3z1jnNSoTF+xHULZtf6HJQMKRNCuKBnn32WlJSUo96/8847SUtLIy0tjTFjxuiQrHmkgLuJzF0lVJXW0XOo81ccO1U9zoqloqiWrD2OmQdYCCFOVVZWFp9//jk33XRTi/eV5URFi+zakINPgJEOLjh07Hg6pkbg5efJrvU5JKSE6h1HCOGCnvj1CXaX7LbrMbuHdmfWGbNOuM3MmTN58sknqaysPOqzF154gXfeeYcBAwbw1FNPERLy1wmzZDlR0WzVFfUc3lZEt0ExbjXHuKfRQLeB0RxMK6SmSjqzCSFcw2effUZkZCT9+/c/6rPp06dz4MAB0tLSiImJ4Z///OdR28hyoqLZ9m3Mx2ZTpAxx3aFjx5MyJJbfVmexb2MBfYa7fuc7IYRznexO2RE2bNjAihUr+OKLL6itraWiooJrr72WRYsWERUVdWS7m2+++S8d3P5MlhMVzbL751wikwIIjfHTO0qLhcf7E57gz56fc/WOIoQQADz++ONkZWVx+PBhFi9ezIgRI1i0aBHAkWfW0NhUfqxn1bKcqGiWoqwqijKrGHplV72jnLLug2JY//E+SnLMhMa635cQIUT7ce+995KWloamaSQnJ/Pqq68etY0sJyqaZc8veXgYNLqcHql3lFPW5fQoNizdz55fchl8cWe94wghxBHDhg1j2LBhR/773XffPek+spyoOCllU+zflE9izzB8/E16xzllvoEmElJC2bexAKVkTLgQov2R5UTbmdwDZVSV1tH19KiTb+ziup4eSWVJLXkHZYETIUT7I8uJtjN7f83H0+RBch/3Gft9PB1SIzAYPdj3a57eUYQQok2QAu6irFYbB7YU0qFvhEsvXNJcJm9PknuHs39LATarTe84Qgjh9qSAu6jsPaXUmi107u++ndf+rnP/SGoqLeTsK9M7ihBCuD0p4C7qwOYCjF4GEnu2nSlIk3qH4WnyYP+WQr2jCCGE25MC7oJsVhsH04pI7hOOp9H9m8+bGE0GknqFcXBrATZZoUwIoaOysjIuu+wyunfvTkpKypHJV0pKSjjvvPPo0qUL5513XrNmUZPlRMUR2fvKGpvP+7Wd5vMmnfo1NqPnSjO6EEJHM2bMYPTo0ezevZtt27YdWVZ03rx5jBw5kn379jFy5EjmzZunc9LjkwLugg5uLcTT5NGmms+bJPUKw2D04GCaNKMLIfRRUVHBunXruPHGG4HGGdSCg4MBWL58OZMnTwZg8uTJfPrpp0ftL8uJimNSNsWhtEISe4bhaWo7zedNTN6eJKSEcnBbIWdd0eWYCwIIIdqPvMceo+53+y4n6pXSneh///u4nx88eJCIiAimTJnCtm3b6N+/P88++yx+fn7k5+cTE9O4cFRMTMxfinMTWU5UHFN+egXm8no6prZ+kL+r6pgaTlVJHUWZVXpHEUK0Qw0NDWzZsoXp06ezdetW/Pz8WtRULsuJimM6lFaEh4dGUq8wvaM4THLvcDQNDqYVEpEYoHccIYSOTnSn7Cjx8fHEx8czcOBAAC677LIjBTwqKorc3FxiYmLIzc0lMvLYfZFkOVFxlEPbContGoy3n1HvKA7jE2AipnMwh7bJc3AhhPNFR0eTkJDAnj17APjuu+/o0aMHABdddBFvv/02AG+//Tbjx48/an9ZTlQcpbywmtK8anoOjdM7isMl9w7nx0/2U1FcQ2CYj95xhBDtzPPPP88111xDfX09HTt25M033wRg9uzZXHHFFSxYsIDExEQ+/vjjo/Zt18uJappmADYB2UqpsXpkcEWHfysGaBNzn59Mcp8wfvxkP+nbi+k9LF7vOEKIdiY1NZVNmzYd9X5YWBjffffdCfdt78uJzgB+1+ncLuvw9iJCYvwIimj7d6Qh0X4ERfpw+LcivaMIIYTDOWI5UaffgWuaFg9cCDwK3OXs87uqupoGcvaWkXpegt5RnCa5Tzjb12ZRX9uAyVue5ggh2q4rr7ySK6+80q7H1OMO/BngXkCWpPqTrN9LsNkUSb3afvN5k+ReYdgaFFm7HdNDUwgh2jJNKefNSa1p2lhgjFLqVk3ThgF3H+sZuKZpU4GpAFFRUf3t1WMPoKqqCn9/f7sdz16yf7FRkQXdL9bQPFx3chN7Xj+bVbFnmSIoCWJPbz8DIlz1Z9BdyPVrHVe4fkFBQXTu3FnXDKfKarViMDhukq39+/dTXl7+l/eGDx++WSl11ENzZ7dbnglcpGnaGMAbCNQ0bZFS6to/b6SUeg14DWDAgAFq2LBhdguwdu1a7Hk8e1BK8faXG+jQO5jhI1o/vZ4j2fv61e79jcLMSs45Z0i7mZXNFX8G3Ylcv9Zxhev3+++/ExDgnnNAVFZWOjS7t7c3p512WrO2deptj1LqX0qpeKVUMnAVsPrvxbs9Ks6uwlxeT1Kvtjf3+ckk9QqjqqSOklyz3lGEEMKttJ92SxeWvqNx+Fhiz7Y7+9rxNP2ZM3aU6JxECNGePP300/Ts2ZNevXpx9dVXU1tbC8hyos2ilForY8AbZe4qISzeH78gL72jOF1AqDehsX5k/l6sdxQhRDuRnZ3Nc889x6ZNm9ixYwdWq/XI7GiynKhotvraBnIPlJPYo/01nzdJSAklZ185lnqr3lGEEO1EQ0MDNTU1NDQ0UF1dTWxsLCDLiYoWyNlXhs2qSGjPBbxHKNu+yyR3X1m7fIwgRHv2w0d77b4yYXiCP0Ov6Hrcz+Pi4rj77rtJTEzEx8eHUaNGMWrUKABZTlQ0X+auEjyNHsR0av28uO4qtkswBk8PMn6X5+BCCMcrLS1l+fLlHDp0iJycHMxmM4sWLWr2/rKcqAAg8/cSYrsG42l03LhCV2c0GYjpHETmLingTZTVirW0FGtFJbbqamxmM6q25i/baEYjHr6+ePj54eHnhyEsDA+v9tePQri3E90pO8q3335Lhw4diIiIAOCSSy7hxx9/5Nprr3Wr5USlgOuoqrSW0rxqepwVq3cU3SWkhPLTsgOYy+vaTWc+m9lM7Z691B86RH16OvWHD2PJyqKhsJCGkhKwtrxPgEdgIJ4RERijozElJWFKTsKUnIxX1654RkW1m7H2QpxIYmIiP//8M9XV1fj4+PDdd98dWVykaTnR2bNnn3A50VdffZVJkyZRUFDAmjVrmDhx4l+WEx08eDAWi4W9e/fSs2fPI8uJDho0SJYTbQuaphCN7x6icxL9NV2DrN2ldBsYrXMa+1NWK7W7d1OzaRNB33zLgSeepP7wYWiaCdHTE1N8PMaEBLxSuuMZEYFnRASGwKDGO2xfXzy8vcDjf0+9VH194915dTXWykqsxcU0FBTSUFiIJTeX8hUrsFX979miITQU75QUvHv1wndAf3xOOw2DzGgm2qGBAwdy2WWX0a9fPzw9PTnttNOYOnUqIMuJimbK2l2KT4CRsFj5JRqeEICXrydZe9pOAa/PyKBq7fdUbVhPzeYtR4qpMSQEU/9+BI4bi3dKCl6dOmGMjUXztO9fR6UU1pIS6g8fpnb3bmp37aJ21+8UL1hA8auvgocH3ikp+A0ehP855+CTmopmNNo1gxCuau7cucydO/eo991pOVEp4DpRSpG1u4S4biEuPfe5s3h4aMR1CyFrdwlKKbds6lVKUbtjBxVfrqJq9erGO2zAlJxM4NgL8R1wOr4D+rNh9256OWEqS03T8AwLwzMsDN/+/Y+8b6uupmbbNqo3bqJ640aK33qb4jcW4BEQgN9ZZxJ4/vn4n3MOHj5tf1lbIZzl888/5/HHH6ehoYGkpCTeeuutVh9TCrhOyvKrMZfXE99Nms+bxHcL4eDWQsoLawiO9NU7TrPVHThA+afLqVi1CktmJhiN+A0cSMg11+B/ztmYEhP/usPu3foE/YOHry9+gwfjN3gwANaqKsw//kjVunVUrf2eyi9Xofn6EjBsGIHjxuI/dKjdWweEaG8csZyo/K3Uyf+ef7ff8d9/l5DSeC2ydpe6fAG3mc1UfPklZUuWUpOWBgYDfoMHEz5tGgHnjsRgh+dbzmLw9ydw1CgCR41CWa1Ub9xIxRdfUvn111R88QWeEREEXXwxwZdegikpSe+4oo1w15Y2R2rp6qBSwHWSvbeUgFBvgiKkmbJJUKQPfkEmcvaW0uvsOL3jHFN9RgYlixZRvvQTbGYzpo4dibznHoLGX4RnuPuv5a4ZDPgNGoTfoEFE3/cfqr7/nrIlSyl+4w2KX3sNv7POInTyJPzOPBPNQ6aREKfG29ub4uJiwsLCpIj/QSlFcXEx3t7ezd5HCrgOlE2RvbeM5F4y69ifaZpGbNcQsveUuty38+rNmyle+CZVq1eDwUDgBRcQcvXV+JyW6lI57UkzGgk491wCzj0XS34B5Z8spfT9D8i8eSqmDh0InTyJoIsvlrHnosXi4+PJysqisLBQ7ygtVltb26Ii2xLe3t7Ex8c3e3sp4DooyTVTW2Uhtqs8//67uK7B7NuYT1l+NSHRfrpmUUpR/fPPFL30MtUbN2IIDiZs2i2EXHU1xqhjT+7QVhmjIgmfPp2wG2+k4quvKHn7HfIenEvRSy8TduMNBF9xhXR6E81mNBrp0KGD3jFOydq1a5u9XrejSQHXQfbexuffcV2DdU7ieuL+6NSXvbdM1wJu/vkXCp95hpq0NDwjI4n6978Ivvzydl+kNJOJoHHjCBw7lupffqHopZfJf3weRa++RthNNxFyzUS5IxfCSeQhlg6y95QREOZNYHj7LgbHEhThg1+wF9l7HDP14MnU7NxJxo03kXH99Vjy84l+8AE6ffsNoZMmtfvi/WeapuE3aBBJ77xN0nuL8O7enYInn+TA+aMpW7oU1dCgd0Qh2jwp4E6mbIqcfWXEdZG772PRNI24rsFk7ytrcY/M1rDk5pJ99z0cvvQyanfuJHLWLDqt+pKQq67Cw2RyWg535Nu/P4kLF5D41lt4RkaSO+c/HBw/gaof1usdTYg2TQq4k5Xkmak1W4iV5vPjiu0STE1FPeUFNSffuJVstbUUvvQSB8ZcSOU33xA27RY6ffM1YVOul6bgFvIbNJDkDxcT99yzqAYLmTffTOb0W6lPT9c7mhBtkhRwJ8vdVwY0FilxbE3XJuePa+UolWvWcPDCsRQ99zz+Q4fS8fPPiZw5E0NAgEPP25ZpmkbgqFF0XLmSyLv/SfUvv3Bg7DgK5j+NrbZW73hCtClSwJ0sZ18ZfkEmef59AsFRvvgEGB1WwC0FBWTNvJOs6bei+XiT+NZbxD/3LKZ41xx77o48TCbCbrqJjqu+JGjMBRS/9hoHLxqP+W/LLAohTp0UcCdSqvH5d2yX4DY7dtgeNE0jtkuw3Qu4UorSjz7i4IVjqVq9moiZM+j4ySf4DRpo1/OI/zFGRhL7xBMkvvUmaJAx5QZyZs3GWubY1hUh2gMp4E5UUVSLubxems+bIbZLMJUltVSW2KfZ1ZKbS+ZNN5N3/wN49+hBxxXLCZ82DU06qDmF36BBdFyxgrDp0yj//HMOjBtH5dq1escSwq1JAXei3P2Ndx0xnaWAn0zTNWq6ZqdKKUXZp59y8KLxVG/dSvSDD5D41puYkpPtkFK0hIeXF5EzZtDh44/wDAkla9p0cubMwVpZqXc0IdySFHAnyt1fhpevJ6Ex+s4w5g7C4vwxehvI3V9+ysewlpeTPWMmubP/hVe3rnT8dBkhV10ljy905p2SQvKSjwm75RbKl33KofETqN6yVe9YQrgdKeBOlHugnOhOQbL+dzN4eGhEdwwi98Cp3YFXb97MwQkXU7l6NZH33E3S228fvayn0I2HyUTknTNJfv898PAg/brrKHrlFZTVqnc0IdyGFHAnqamqpzSvmphO7rPMpN5iOgVRnGOmrtrS7H2U1UrhSy+Rft0kNKOR5A/eJ+zGG9EMBgcmFafKJzWVDss+IXD0aAqfeZaMKTdgyS/QO5YQbkEKuJPkHWhsCo7pJM+/myumUxAoyDtY0aztG0pLybxlGkXPPU/ghRfS4ZOl+PTu7eCUorUMAQHE/t9/iXnsMWq2b+fQpZdi/vVXvWMJ4fKkgDtJ7oFyPDw1IpNlkpDmiurQ+Lghpxkd2Wq27+DwpZdR/csvRD80l9gnn8Dg7++ElMIeNE0j+JKL6fDRhxgCAsiYcgPFCxY6dTpdIdyNFHAnyTtQTkRCAJ5GacptLqOXgYgE/yOtF8dTtmQJ6RMnApD0/vuEXHGFdFRzU15dupD88UcEnHsuBf/9L9kzZmIzm/WOJYRLkgLuBNYGGwXplUTL8+8Wi+4YREF6BVar7ajPVEMDeY8+Ru5/7sP3jDNIXroEn969dEgp7Mng70/cM08TOWsWld9+y+GJ11Cfla13LCFcjhRwJyjMrMTaYCOmoxTwloruGERDvY3irKq/vG8tLydz6i2UvvsuoZMnk/DqK3iGhOiUUtibpmmETbmehNdfw5Kby+HLL6d60ya9YwnhUqSAO0H+H52woqWAt1hTq8WfO7LVHTrE4SuuxLxxIzGPPkLUv2ajeXrqFVE4kP+ZZ5L84WIMQUGkT7mBsqVL9Y4khMuQAu4EuQfKCQj1xi9YlqdsKf8QL/yCvcg72PgcvHrLFtKvnoi1ooKkt94k+NJLdU4oHM2rQweSP/oQv9NPJ3fOf/BbsUI6twmBFHCnyD9UTnTHQL1juCVN04juGEjewXIqVq0i4/opGIKCSF78Ab79++sdTziJITCQhFdfIejSS/D/4ktyZ89G1dfrHUsIXUm7o4NVldZSVVpHlDSfn7LojkEc2FLIwXsfILhnT+JfelGed7dDmtFIzCOPkG2xwPIVWPLyiX/heVm/XbRbcgfuYE3PbqM7SAE/FUopTD9/AYDl7AkkvrlQinc7pmka5jFjiH1iHtWbN5M+eTINxcV6xxJCF1LAHSz/cAUGTw/CE2RSkZZSVit599+P7f2X0LBhGXEZHt7eescSLiBo/HgSXn6J+oOHSJdhZqKdkgLuYPmHyolI9MfgKZe6JWz19WTPvJOyj5cQecuNhCcFUXBYlp0U/+M/dCiJCxfSUFpK+rH+m9MAACAASURBVMSJ1O3bp3ckIZxKqooDWa02CtMriUqW5vOWsNXWknXrbVR+8w1R/5pN5IwZRHcIIj+9EptNeh+L//HtdxpJ774LSpF+3SRqd+3SO5IQTiMF3IFKss00WGxESQ/0ZrNVV5N5yzTMGzYQ8+gjhE6eDEBUh0Aa6qyU5Mi0muKvvLt1Jem9RWi+PqRfP4Wa337TO5IQTiEF3IGaxi5HJUsBbw5rVRUZN0+leuNGYp984i9jvKM6NF7D/EMnnhddtE+mxESS330XQ1AQGVNuoHrLFr0jCeFwUsAdqCC9Ap8AIwFh0vHqZKyVlWTceCM127YRN/8pgsaN+8vnQRE+ePl5UnC4eUuLivbHGBdH0qJ38YyIIOPGm2RJUtHmSQF3oPzDlUQlB8rKWCdhrTKTefNUanfuIv7ZZwgcPfqobTRNIyopkHzpyCZOwBgVRdK772CMjSVz2nSqN2/WO5IQDiMF3EHqaxoozTMTKc3nJ2Qzm8m85RZqtm8nbv5TBIwcedxtI5MDKcmpwlJndWJC4W48IyJIfHMhxshIMm+eSk1amt6RhHAIKeAOUpBRCQop4Cdgq64mc9p0atLSiHvq/wgcNeqE20clB6IUFGbIXbg4MWNkJIlvv4UhIpyMm26mZvt2vSMJYXdSwB2k6VltVJIU8GOx1deT9Y87qN68mdgnnjhms/nfNX0ZKkiX5+Di5IxRUSS99RaG4GAybryJ2j179I4khF1JAXeQgsMVBEb44O1v1DuKy1ENDeTcfU/jULGHHyZo7IXN2s830IR/qBf50pFNNJMxJobEt97Cw8eHjBtvoj49Xe9IQtiNFHAHyT9cQVSSLLLwd8pmI/f+B6j8+mui/v0vgi+9pEX7RyUHSk900SKm+DgSFy4Aq5WMKTdgycvTO5IQduHUAq5pmremab9qmrZN07SdmqbNdeb5naW6op6q0jp5/v03SikKnniS8k8+Ifz22wmdNKnFx4hMCqSiqJZas8UBCUVb5dWpEwlvvI61ooKMG26koaRE70hCtJqz78DrgBFKqb5AKjBa07RBTs7gcE3PaCPlDvwvil9/g5K33yZk0nWE33brKR2j6ZoWpktHNtEyPj17kvDKy1iys8m8ZRq26mq9IwnRKk4t4KpR1R//afzj1eYmty7MqAQNwhOkgDcpW/YphfPnEzh2LFGzZ5/y2PiIxMZrWpAhzeii5XwHDCDu6fnU7txJ1p13oizSkiPcl6aUc+unpmkGYDPQGXhRKTXrGNtMBaYCREVF9V+8eLHdzl9VVYW/v2OX9sxYZ6O+CjqPaXtdDE7l+pl27iT4xZeo79qVsttvA0/PVmXY97kNryBIPMs9r68zfgbbMntcP58f1hP43nvUDB5MxaTroB1NtiQ/f62jx/UbPnz4ZqXUgL+/37rfpKdAKWUFUjVNCwaWaZrWSym142/bvAa8BjBgwAA1bNgwu51/7dq12PN4x/Lml+tJTgll2LAeDj2PHlp6/Wq2byd9wUJM3brR7d13MNjhB7/+wE5yD5QxbNiZrT6WHpzxM9iW2eX6DRtGYXAwRS++SHyfPkTeOdMu2dyB/Py1jitdP91uYZRSZcBa4OQDgN2IuayO6vJ6IuT5N/VZ2WROm45nSAiJr71ql+INjc/Bq0rqqK6ot8vxRPsUfvttBF9+OcWvvkrpRx/pHUeIFnN2L/SIP+680TTNBzgX2O3MDI52pANbYvsu4NaKCjKn3YKyWEh47VU8IyLsduymjmwyoYtoDU3TiH7gfvyGDiVv7kNUbdigdyQhWsTZd+AxwBpN034DNgLfKKU+c3IGhyrIqERr5x3YlMVC1owZ1KdnEP/cc3h16mTX4zdd26JM6YkuWkfz9CTu6fl4depE9oyZ1O7dq3ckIZrN2b3Qf1NKnaaU6qOU6qWUesiZ53eGooxKgqP9MHoZ9I6iC6UUuXPnUv3Tz8Q89BB+gwba/Rwmb0+Co3wpzKg6+cZCnITB35+EV1/Bw8eHzGnTaCgs1DuSEM3int14XVhhRiURie23h2fJwjcpX7KU8FunE3zxBIedJyLBXxY1EXZjjIkh/pWXsZaWkXnb7djq6vSOJMRJSQG3I3N5HebyeiIT2+cMbFXff0/B//0fAaNHE/6Pfzj0XOGJAVSW1FJbJeN4hX349OxJ3H+fpPa338i97z6cPcRWiJaSAm5HTXeE7fEOvO7AAbL/eTdeKd2JffyxU56opbmaOgnKXbiwp4BzzyVixh1UrFhJyYIFescR4oSkgNtRUzFpbx3YrGVlZE6/Fc3bm4QXX8TDx8fh52y6xjIjm7C3sGnTCBxzAQVPzadyzRq94whxXFLA7agwo5KgSB9M3k6fH0c3qqGB7LvuoiE3l/jnn8MYE+OU83r7GQkM95aObMLuNE0j5tFH8e7Rg5y776HuwAG9IwlxTMetNJqmvXOKx7xPKdUuF90tyqwiqmP7ev5dMP9pzD/+RMxjj+F72mlOPXd4QoAMJRMO4eHjQ/yLL3Do0svIuu12kj/+CENA+2pZE67vRLeK1wJ5NK4g1lyJwDNAuyvgtWYLlSW19DonTu8oTlPxxReULFxIyMSrCb7kYqefPyLBn4NbC6mvacDk035aPYRzGKOjiX/madKn3EDOrNnEv/A8moc0WgrXcbKfxglKqQ7NeQFdgPazIsDfNN0Jhie0jw5stXv2kjPnP/j060fU7Nm6ZDgyoUuWNKMLx/A9/XSiZs2iavVqil55Re84QvzFiQr4d0B5C45l/WOfdtmrqDCzsYiEx7f9ZjZreTlZ//gHBn9/4p99Bs1k0iVHxB8FvFCa0YUDhVx7DUHjL6Lo+ReoXLtW7zhCHHHcAq6UOk8ptae5B/pjre/zlFL77RPNvRRlVuIXZMI3UJ9i5jQ2GzmzZmPJzSXuuWftOsd5S/kGmfAJMMpzcOFQmqYRPXcuXindybl3FvVZWXpHEgKQXuh2U5RVRXg7WMDE95tvqFq7lqh773V6p7W/0zSNiIQAaUIXDufh7U38s8+CUmTPmCkztQmX0OwCrmmap6ZpQzVNm6hp2g1/fzkypKtrqLdSmld9pEm3rareuBH/5SsIuGA0Iddeo3ccoPE5eEmOGavFpncU0caZEhKIfWIetTt3kj9vnt5xhDhhL/QjNE3rBywD4jl2RzUFLLRjLrdSnGNG2RTh8W23A1tDYSFZd92FNSKcmIcfcfhMa80VnuCPzaooyTO3+S9QQn8BI0YQeuMNlCxYiG+//gSNG6t3JNGONfcO/BWgCpgAdAM6/O3V0SHp3ERb74GurFay774HW2UV5VOnYvD30zvSEU1fmooypRldOEfkzJn49O9P7gMPyCQvQlfNLeA9gHuUUiuVUvuUUul/fzkypKsrzqrC6GUgMMzxU4jqoeiVV6j+5Rei77+fhjjXGuceFOmLp9GDYnkOLpxEMxqJmz8fD29vsu+8C1ttrd6RRDvV3AK+F3Cd2y4XU5RdRXi8P5qHazQr25P5118pevElgsZfpMtkLSfj4aERGudPUbb0RBfOY4yKJPaJedTt3SvPw4VumlvA/w38R9O0REeGcUfKpijKqiKsDT7/bigtJeeeezElJBB13/16xzmu8AR/ijKrZPlH4VT+Q4cSeuMNlC3+kIpVX+kdR7RDzerEppRapWnaMGCfpml7gdKjN1Hn2DucO6gorsVSa21zHdiUUuT+699YS0pI+HCxSz33/rvwOH92/ZBDVWkdAaHeescR7UjkzJlUb9pE7n334d2rJ6b4eL0jiXakub3QZwP3AoU0zrRmdWQod1KU9UcHtjY2A1vpO+80jveeMwfvHj30jnNCRzqyZVVJAW8Gi9VCeX05VfVVVFkaX2aLmXprfePLVo/FajlqP08PT0wGEyYPEyaDCV9PX/xN/vgb/fEz+hHkFYS3Z/u6/prRSNxTT3Ho4kvI/uc/SV60CM1o1DuWaCeauwLETOBV4HallBTvPynKqkLTIDTOde9QW6p2zx4K/u8p/IcPd5nx3ifS9PiiOKuSDn3CdU6jr/K6cjIrM8muyibPnEd+dT555jyKaoooqS2hpLaEynrH9Rfw9fQlxDuEMO8wwnzCiPaLJso3imi/aGL8YkgMTCTMO8xlhiHagyk+npiH5pJ9510UvfwyEXfcoXck0U40t4D7Ah9L8T5acVYVQZG+GE0GvaPYha22lpy778YjOIiYR11nvPeJmLw9CQz3pijLrHcUp7ApG9mV2ewv28+B8gPsL9tPenk6mVWZlNf9dfkCH08fonyjiPSNJCU0hVDvUEK8Qwj2Csbf5E+AMQA/ox9+Rj+8DF4YDUZMHiaMBiPan6Z8UCgabA1H7tDrrfXUNNRQWV+J2WKmsr6SivoKimuKKa0rpaSmhMzKTDblbaLSUnlUpsSARBIDE+kU3IlOwZ3oHNSZpKAkjB7uefcaeMEFVK37gaJXXsXvzDPx7d9f70iiHWhuAf8SGAysdmAWt1ScXUVEYttZA7zg/56ibt9+El5/Hc/QUL3jNFtYnD/F2W1vKJnVZmV/2X52Fe/i95Lf+b34d/aU7qGmoebINtF+0SQHJnN+0vkkBiaSEJBAnH8c0X7RBJoCdf8SZraYyTfnk12VTWZlJpmVmWRUZrCnZA/fpn+LorHzoaeHJ12Cu9AjrAcpoSmkhKXQPbQ7JoN7rC8QNWcO1Zs2kXPPvXRY/qmsHy4crrkF/BngrT9+Eazi6E5sKKUO2jGXW6ivbaCiqJaUIbF6R7GLqu+/p3TRIkInT8J/6Fl6x2mR8Hh/Dv9WhKXe6tatIZX1laQVpJFWmMa2gm1sL9pOdUM10Ng83T20O5d2uZSuIV3pGNyRTkGd8De5dgdKP6MfHYM70jH46PmeahtqOVxxmP1l+9lXuo/dJbv5LuM7lu5bCoDRw0hKWAqpEan0jehLv6h+zo7fbAZ/P+L++ySHr7mWvIceJu6/T+odSbRxzS3gG/7458PAQ8fZxn1/a56ikpzGJtu2MISsoaSEnH/PwatrVyLuukvvOC0WFu+PUlCaayYyyX1aRKot1WzO38zGvI1szNvIrpJd2JQNg2aga0hXxnUaR9+IvvQK70VSYBIeWttaf8jb05vuod3pHtr9yHtKKfLMeews3slvhb+RVpjG4t2LeWfXOwBEG6NZ//N6zog+gzOizyDYO1iv+EfxSU0l/LZbKXruefzPPlumWhUO1dwCfgMgg2z/pmkVrDA378CmlCLvgQewVVQQu3AhHl5eekdqsbC4//VEd+UCrpQiuz6bN3e8yYbsDWwp2ILFZsHTw5M+4X24uffNDIgeQJ/wPvgaffWOqwtN04jxjyHGP4Zzk84FGnvO7yrZxaa8TXy16ytWHFjBh3s+REOjV3gvhsQO4ay4s+gV3gtPj+b+WnOM8KlTMf+wnryHH8b39AEYo6N1zSParuaOA3/LwTncUnF2FSZvg9sPXSpfvpzKb74l8p578O7WVe84pyQo3AdPk4dLPge3WC1szN/Imow1rM1aS545D3Khc3Bnrkm5hsGxgzkt8jR8PNvmVLz2YDQY6RvRl74RfelU3Ikzzz6TnUU7+Sn3JzZkb+D17a/z6m+vEuQVxDnx5zA8YThDYofo8iVI8/Qkdt7jHJxwMbn/nkPCgjd074cg2iZ9v6q6ueLsxhnY3PkvpyUnh/xHHsVnQH9Cr5+sd5xTpnloLtWRrd5az4bsDXyV/hXrMtdRaanE2+DN4NjBjPAawZQRU4j2kzuzU2X0MJIamUpqZCrT+06nvK6cn3J+4vus71mbuZYVB1Zg8jAxOHYw5yefz7CEYQSYnNepzJSURNSse8l7cC6lH3xA6MSJTju3aD+OW8A1TVsIPKyUOtScA2mNVWwB8KBSKsNO+VyWUorirCq6nuG+v4SVzUbOv+eAzUbs44+jGdy7G0NYnD8HthaglNLlS5XFZuGnnJ9YdWgVazLXUGWpIsgriHOTzmV4wnAGxQ7Cx9OHtWvXSvG2syCvIEZ3GM3oDqNpsDWwtWArqzNW823Gt3yf9T1GDyNnxp3J6OTRDE8Y7pQ78+Arr6Ty2+8oePK/+A8Zgik52eHnFO3Lie7AJwMvA80q4DTOqz4ZeAFo8wW8qrSO+lqrW3dgK130HtU//0z0ww9hSkjQO06rhcX5sWt9A9UV9fgFOec5vlKKXcW7WHlwJV8e+pKS2hICTYGcl3Qe5yefzxkxZ7jt2GZ35enhyenRp3N69Oncc/o9bC/azleHv+Lrw1+zNnMtvp6+nJt0LmM7juWM6DMweDjmi6umacQ8+ggHLxpPzqzZJL23CM1TGj2F/Zzop0kDXtY0raKZx3LfduRT0NRUGxrrnh3Y6g8fpmD+fPzOOZvgyy7TO45dhMX+MSNbdpXDC3hJbQkrD6xk2b5lHCg/gMnDxDkJ5zCu4zjOijsLo0GKtivw0DyOPDu/e8DdbMnfwmcHP+Prw1+z4sAKIn0jmdB5Ahd3vpj4APvPY26MiiL6vvvIuftuSt56i7CbbrL7OUT7daICvo7GnuctKczrgHaxrmNTAQ9zwwKubDZy5vwHzWQi5qGH3foZ/p81TWdbnG0msUeY3Y9vUzZ+yf2FpfuW8l3GdzTYGugb0ZcHBj/AqORRBJpct/e7aCzmA6IHMCB6AP8a+C/WZq7l0/2f8sb2N3j9t9cZFDOIS7pewsiEkXb9AhZ44Rgqv/qKwueex3/4cLw6dbLbsUX7dtwCrpQa5sQcbqc424x/iBdevu53p1W6aBE1mzcT8/jjGKMi9Y5jNz7+JnwDTZTYuSNbZX0lKw6sYPHuxRyuOEyQVxBXdbuKS7tcSueQznY9l3AOL4MX5yefz/nJ55NnzmPZ/mUs27eMe76/h3CfcC7vejmXdb2MSN/W//3QNI3oB+7n4IVjyfn3v0l+/323728iXIM8kDlFJTnmI2OP3Ul9ejoF85/G75yzCZowXu84dhcW50dxjn3mRD9YfpD3dr3HyoMrqWmooU9EHx476zHOTz7fbab3FCcX7RfN9L7Tmdp7Kj/m/MjiPYt5ZdsrvP7b64xIHMF1Pa4jNTK1VefwDA8n6s9N6TfeaKf0oj2TAn4KrFYbpXlmEnu6z1zh0NR0PgfNaCTmoYfaTNP5n4XG+bPj+2xsNoWHR8v/fEopfs37lXd2vcO6rHWYPEyM6TiGq7pfRc+wng5ILFyFwcPA0PihDI0fSmZlJh/t+YhP9n3C1+lf0yeiD5N6TGJk4shTniimsSl9FYXPPof/sGHSlC5aTQr4KSjLr8ZmVW53B176/gfUbNpMzKOPYIyK0juOQ4TF+mO12KgorCE4qvlDhaw2K1+nf83CHQvZXbKbUO9Qbk29lSu7XUmot3t9UROtlxCQwD8H/JPpfaez/MBy3t31Lnd/fzdx/nFc1+M6LulySYsn3tE0jej77+fg2HHkzvlPY690aUoXrdC2JlZ2kqY50N2pB7olJ4fC+fPxO/NMgi65RO84DhN2pCNb856D11vrWbp3KRd9ehH3rruXOmsdDw5+kK8v+5rpfadL8W7nfI2+XN39alZOWMkzw58h0jeSeb/OY/TS0byx/Y0Wr63uGRFB1L//RU1aGqXvf+Cg1KK9kDvwU1CSY0bz0AiJdo+5qpVS5D74IEopoufObZNN501CYvxAg5JcMydqoKyz1rFk7xIW7lhIQXUBPcJ68MywZxieOLzNLRgiWs/gYWBk4khGJo5kc/5mXt/+Os9ueZYF2xcwMWUik3pMIsgrqFnHCrzoIspXfkbB008TMGI4xrg4B6cXbdVxf1NpmvYPZwZxJ8XZVQRF+OBpdI/mr4rPPsO87gciZ87AFN+2f1kYTQYCw30ozj52R7Z6az0f7P6AMZ+MYd6v80gISODV815l8YWLGZk0Uoq3OKn+Uf155dxX+HDshwyOHcxrv73G6KWjeTHtRSrqTz5thqZpxMx9EIDcB+eilKwTJU7NiX5bPaNp2jpN02SczN+U5JrdZvx3Q0kJ+Y8+hnffPoRce63ecZwiLNaPkty/FnCLzcJHez5izCdjeOyXx4j3j+eNUW/w5vlvMiR2SJtulRCO0SOsB/OHzWfJuCUMihnEK9teYfSS0by87WXMlhOPhDDGxRE5cybmH36gYuVKJyUWbc2JCvhIIAbYpmna3Zr8hgOgod5KeWGN2zz/zn98HlazmdhHHmk3HWZCY/woz6/GarFhUzZWHVrFhE8n8PDPDxPjF8Pro17nrdFvMTBmoBRu0WrdQrvx9PCn+Xjcx5wefTovpb3EmE/G8N7v71FvrT/ufiHXTMSnb1/yH3uchpISJyYWbcVxC7hSai3Qm8a5zR8DftI0rYeTcrms0rxqUBAa6/o90KvWb6Bi5UrCb74Zry5d9I7jNKFxfthsirXbf+Sqz67innX3YDKYeGHEC7xzwTsMihkkhVvYXffQ7jw74lneH/M+nYM7M+/XeYxbNo4VB1ZgU7ajttcMBmIeeRir2UzBE0/qkFi4uxM+8FNK1SqlZgGDABOwRdO0xZqmvfO319tOSesCSnLcYw50W20teXPnYkpOJuyWqXrHcSqzX+PdzPNrXqOivoLHznqMJeOWcE7COVK4hcP1jujNG6Pe4NXzGtcnn7N+Dld9dhUb8zYeta1Xly6E3XgD5cuXY/75Zx3SCnfW3B47+4E0Gov40OO82oXiHDMeBo2gyJaNAXW2opdfwZKZSfSDD+Lh5ZyVufRWXFPMwz89zA2/XosNK2OCLmbFhBWM6zTOYStOCXEsmqYxJHYIi8cu5vGhj1NaV8oNX93AHavv4HD54b9sGz5tGsbERPIeeBBbXZ0+gYVbOmkB1zTtImAXMAG4WSkVp5Tq8LdXR4cndREluWZCon0xGFy3t3Ldvn0UL1hA0IQJ+A0aqHcch7PYLLyz8x3GLhvLJ/s+4YqUywmO9CXe0kmmPBW68tA8GNtxLCsnrGRGvxn8kvsLFy+/mKc2PUVVfWNrnoe3N9EP3E99ejrFr76mc2LhTk40jCxC07QPgWXAFqCXUmqB05K5qJIcc+NYYxelbDZyH3gQg78/kbPu1TuOw/2Y8yOXrbiM/276L6mRqXwy/hP+NfBfRMQFHtUTXQi9eHt6c1Pvm/j8ks+5qPNFvL3zbcZ9+r/n4/5nnknguHEUvf46dQcP6h1XuIkT3UbuBkYAk5RSFymlcpyUyWVZ6qxUFtcS6sIFvGzJEmq2bCHy3nvxDAnRO47D5FblMnPNTG755hYsNgsvjHiBl0a+RIegDkDjhC4VhTU0WKw6JxXif8J9wpk7ZC7vX/g+sX6xzFk/h0lfTmJ3yW6iZs/Cw9eXvPsfkLHhollOVMBXAz2VUu85K4yrK81z7SlUG0pLKXxqPr4DBhB08QS94ziExWbhrR1vMX75eDZkb2BGvxl8Ov7Tozqohcb6oVTjvPVCuJpe4b14d8y7PHLmI2RWZnLlZ1cy/+ACgmfcRvWmTVSsWKF3ROEGTrQe+OX2PpmmaQnAO0A0YANeU0o9a+/zOEpTk6yr3oEXPPUUVrOZqPvva5O9rdMK0njo54fYV7qPYfHDmD1wNnH+x55Zrun/UUmumfD4AGfGFKJZPDQPxncez7CEYTy35TkW7VrE197h/F+3JPKf/C/+w4djCAzUO6ZwYc7uidUA/FMplULj0LTb3GlseUlTD/QI1+uBXr11K+VLlhI6aRLeXbvqHceuquqreOTnR7juy+uorK/k2eHP8vzI549bvAGCI33RPLQjC88I4aqCvIK4b/B9LBqziBDfMO4fkoWlpJj0p+bpHU24OKcWcKVUrlJqyx//Xgn8DrjN5NyluWaCo3zxcLEe6KqhgbyHHsYzKoqI227VO45dfZ/5PROWT+CjPR9xbcq1LB+/nBGJI066n8HoQVCED6W50oQu3EOfiD58cOEHXDr2Hr7rb6Tmo2Ws+OKZY04CIwSApldnCU3TkoF1NPZur/jbZ1OBqQBRUVH9Fy9ebLfzVlVV4e9/arOo7V1pwycMEoa4VgH3WbOGwA8/ouzmm6nr38+h52rN9WuJSmslS0qWsKV6CzHGGCaGTSTZK7lFx8hYb6OuHLpc6Fr/v5x1Dduq9nD9SsozSXjoCXKCrLx1U1eujriGSGOkXY7dHq6fI+lx/YYPH75ZKTXgqA+UUk5/Af7AZuCSk23bv39/ZU9r1qw5pf3qaxvUC9O+U79+dtCueVrLUliodvcfoNJvuFHZbDaHn+9Ur19LfHP4G3X24rNV6jup6qW0l1R9Q/0pHefn5QfUi9O+Uw31VjsnbB1nXMO2rL1cv9Jln6pd3bqrWf/spwa8O0C9u/NdZbW1/me5vVw/R9Hj+gGb1DHqo9NvTTRNMwJLgfeUUp84+/ynqiy/cQ70kGjX6sBW8NR8bHV1RP1njtt3XCurLePedfdy59o7ifKN4qOxHzG973SMBuMpHS80prEneqn0RBduKGj8Rfj078/1642cHdSPJzY+wZRVU8isyNQ7mnARTi3gf6xotgD4XSk135nnbi1X7IFek5ZG+bJlhE2ehFeHDnrHaZWmZ93fpH/D7am3896F79ElpHULsITE+AL/G/4nhDvRNI3o/8xBVVRyz7YEHjnzEfaV7uPSlZeyePdiGSsunH4HfiZwHTBC07S0P15jnJzhlJTmmtE8XGcOdGWzkffIo3hGRBA2bbrecU5ZtaWauT/N5fbVtxPmE8biCxdzS99bMHqc2l33nwVH+qJpjf/vhHBH3ikpBF95BaUfLOZ8awrLxi+jX1Q/Hv3lUaZ/N53C6kK9IwodObsX+nqllKaU6qOUSv3j9YUzM5yq0rxqgiJ8MHi6Roeo8k8+oXbHDiLvvQeDv+u0CrTEb4W/cfnKy1m6dylTek3hgws/oFtoN7sd39NkICDcp3EJWCHcVMQdd2Dw9yf/kUeI9I3k5ZEvM2fgHDblbeKSFZfwXfp3ekcUkbobtwAAIABJREFUOnGNauQGSvPMLtN8bq2ooGD+0/j070/g2LF6x2kxq83KK9teYdKXk7DYLCw4fwF39b/LIQuPhEb7ypzowq15hoQQcedMqjdupHLVKjRN46ruV/HRuI+I8Yth5tqZ3L/hfqot8kW1vZEC3gzWBhtlBTWERPvqHQWAwhdewFpWRrQbdlzLM+dx09c38WLai4xKHsWSi5ZwevTpDjtfSIwfZQXV2Kwylla4r+DLL8erRwr5TzyJrbqxUHcM6sh7Y97jpt438en+T7nysyvZXbJb56TCmaSAN0N5QQ3KplxiFbK6Awcofe99gi+/HO+UFL3jtMjqjNVctvIydhbv5JEzH+GJoU8QaHLsVJEh0X7YGhQVRbUOPY8QjqQZDETPmUNDXh7FC9888r7RYGRGvxm8Pup1zBYzEz+fyHu/vycd3NoJKeDNcGQRExco4PlPPIGHry8RM+7QO0qz1VvrefyXx5mxZgaxfrF8NPYjxnce75TWg6ae6NKMLtydb//+BFwwmuI33sCSm/uXzwbGDGTpRUsZEjuEeb/O447Vd1BeV65TUuEsUsCboamAB0fp24RetW4d5nU/EH7rrXiGhuqapbmyKrOY9OUk3t/9PtemXMuiMYtIDkp22vlD/xi3L0PJRFsQdffdYLNRMP/poz4L8Q7h+RHPM/uM2azPWc/lKy/nt8LfdEgpnEUKeDOU5FbjH+qF0cugWwZlsZA/7wlMSUmEXjNRtxwtsTpjNVd8dgUZFRk8M/wZZp0xyyEd1U7E5OOJX7CXzIku2gRjXByhN0yhYuVKatLSjvpc0zSuSbmGdy94Fw/Ng8mrJrNo1yJpUm+jpIA3Q1l+9ZE7Ob2ULv6Q+oMHiZw1C83k3CLYUg22Bp7a9BQz/r+9+w6PqvgaOP6dLUk2m2x6AqFLVUA60jsIIvhaUZFiAQRBEBR7/2FDpEgTpBdBsYEISC+CSlNpIk2lpve+u/P+cRNAJECSLdnsfJ5nn5DNZu7hJtlzZ+7Mmc0jqRRYieW9ltO5cme3xRNSzl9VY1PKjPBBgzBERHDh3XeR9qtPzqwXXo/ldy6nTYU2vL/7fUZvGU16brqLI1WcTSXw65B2SdKFDLeWULUlJxM3dSrmVq0I6NjBbXHciPiseAb9MIj5h+bTp3YfFvVYRKXASm6NKaScmaQLGaoXopQJOrOZiDGjyf7td1K/+67Q1wX5BjGl4xSebfosm09v5qHVD3Ei+YQLI1WcTSXw60hPzsGaayfYjUvI4mfMwJ6WRuQLz5fqZWO/xv5Kn1V9OBh/kHfavMMrLV5x+ZD51YSU8ycv20ZmSq67Q1EUhwjq3Ru/evWI/Wgi9uzCV1gIIRhQdwCfdvuUtNw0Hlr9EPsy9rkwUsWZVAK/joLJT+5aA577998kLv2M4Pvuw69WLbfEcD1SSpb9sYxH1z2Kj96HxXcsplf1Xu4O66KCiy81kU0pK4ROR9TzY7FeuEDi/AXXfX3Tck1ZfudyaoXUYl78PMbvHo/VbnVBpIozqQR+HQWTn9w1hB774QR0RiMRI4a75fjXk2vL5Y1dbzDu53G0KN+CZXcuc2g5VEe4NBNd3QdXyg7/Zs0I6NKZhFmzsMbHX/f1UeYo5t0+j7aBbVl4eCFPbniS5OxkF0SqOItK4NeRFJOJr78BU2DJN9coqsw9e0hbv56wQU9giIhw+fGvJz4rnsfWPcZXx75iUP1BTO00lSDfIHeH9R/+QT4Y/fQqgStlTuSYMdhzc4n7eOoNvd6oN/JA6AO81eot9sXs46HVD3Es6ZiTo1ScRSXw60i+kEFIOX+X33uWdjsx73+AISqK0IEDXXrsG3Ew/iB9vuvDn0l/Mr79eJ5u/DR6nfuW2V2LEIKQKH81hK6UOb7VqhHy0EMkf/EFOceP3/D33V3zbubePpdsWzZ9v++rNkTxUCqBX0fShUy3DJ+nfr+G7AMHiBg1Cp2pdGxhWmDNqTUMXDsQgzCwsMdCulft7u6QriukvJlktZRMKYPChw1FZzYTM358kb6vYWRDlvVcRvWg6ozaMopPfvtErdTwMCqBX0NOZh6Zqbkur8Bmz80lbuJEfG++maC7erv02NcipWR18mrGbhtL3bC6fHbnZ9QJrePusG5ISDl/0pNyyM1WE3eUssUQEkL4k0+SsXUbGbt2Fel7o8xRzO8xn5439WTqr1N5cceL5NhynBSp4mgqgV9DQfEPV89AT1q6lLyzZ4l67lmErnT8iLKt2YzdNpa1KWvpXb03s7vNJtTPM8q5wqUyuKoXrpRFIY/0xRBdntgPJxRa3KUwvnpf3m3zLiMajWD1ydU8vu5x4rOuPylOcb/SkR1KqYI3e1f2wG2pqSTMmIm5dWvMrVq57LjXUjBZbd1f67gr+C7+1/p/pWJ9d1GoBK6UZTpfXyJHjiT70CFS16wp8vcLIRh862AmtJ/A0cSjPLz6YTW5zQOoBH4NSRcy0ekElgjX3YNOmP0ptpQUIp8d47JjXsuJ5BP0Xd2X48nHmdhxIl2CupTqYjKFCY7wRwi1lEwpuyy9euFbpw5xEydhzy1e0aJuVbsxv8d8rHYr/df0Z9e5og3JK66lEvg1JMdkYokwode75jTlXbhA4sKFWHr1KhV7ff9y/hf6fd+PHFsO826f59Z65iWlN+oIDDepHrhSZgmdjsgxo8k7c4bkZcuK3U7dsLos7bmU8gHlGbZhGF8f+9qBUSqOpBL4NSTHZLr0/nfcxx+D3U7EyJEuO2ZhVp5YyZANQ4j0j2Rpz6XUDa/r7pBKLCRKbWqilG3mNm3wb9GC+OkzsKWlFbudcuZyLOy+kGblmvHaztf4eP/HaoZ6KaQSeCHsdklybKbL7n/nHDtGytffEPLww/hUrOCSY16NlJJZv8/i5R0v0ySyCQvvWEh0QLTb4nGk4Ch/UmIykXb1RqSUTUIIIp99FltyMgmfzilRWwE+AUzrMo17at5z8T0hz57noEgVR1AJvBBpCVnYrdJlCTx20mR0/v6EPTnEJce7Gpvdxrifx/Hx/o+586Y7mdFlBhYfi9vicbTgKH+seXbSkgrf/EFRPJ2pXl0sd9xB4sKFWOPiStSWUWfkjZZvMLzhcFadXMWIjSPIzFOjWKWFSuCFKJjs5IoiLlm//kr6xo2EPf4YhpAQpx/varKt2Ty79VmWH13Oo/UeZVybcRj1ri8f60wFt0PUfXClrIt4egQyN5f4GTNL3JYQgiENhvBmqzf56fxPPLruUbXMrJRQCbwQl5aQOXcGupSS2I8mog8LI7R/f6ceqzApOSkMWT+Ejf9s5PlmzzO6yWh0ouz9ahSMpqiZ6EpZ51O1KsH33kvSF1+Qe/q0Q9q8p+Y9TO44mZPJJ+m/pj//pP7jkHaV4it779IOkhyTia/ZgCnAueudM37cSeYvvxD+5JPozK4v2RqbGcvAtQP5Pf53Pmj3AY/c8ojLY3AVf4u2qUmK6oErXiD8qWEInU6bHOsg7Su159Pbtb3F+6/pzx+JfzisbaXoVAIvRHJsJiFOvv8t7XbiPvoIY3Q0wX0ecOqxruaf1H/ov6Y/Z9PPMr3zdLpXK/01zUuiYFOT5FiVwJWyzxgVRcgjfUld9R3ZR/90WLsNIhqwoMcCDDoDj659lL0xex3WtlI0KoEXIjkmi+BI5ybwtB9+IPvwYcJHjEDn49rKZkcTj9J/TX8y8jKY020OLaNbuvT47hIU6U9yTJa7w1AUlwh74gl0ZjNxkyc7tN2bgm5iUY9FhJvCGbJ+CFtPb3Vo+8qNUQn8KnKzrWQk5xDkxB64tFqJmzwFnxrVCerdy2nHuZr9sft5dO2jGHQGFnRfQP2I+i49vjuFlPMnLSkba67N3aEoitMZQkIIe+Jx0jdtInP/foe2XT6gPAt6LKBGcA1Gbh7JqhOrHNq+cn0qgV9FSqzWQ3PmEHrKylXknjpFxNNPI/Su20d757mdDP5hMGGmMBb1WMRNwTe57NilQXCkP0hIiVO9cMU7hPbrhz40lLgpUxzftl8oc26fQ5OoJry04yU+P/q5w4+hFE4l8KsouEfqrDXgMjeX+GnT8LvlFgK7dnXKMa5m8z+bGb5xOJUtlZnXfR7lA8q77NilhdrURPE2OrOZsMGDyNz1Exk//ezw9s1GM9O7TKddxXa8/dPbLDi0wOHHUK5OJfCrSI7JBAFBTtrEJPmrr8g7e5aIUSNdtjHImlNreGbLM9QJrcPc2+cSbgp3yXFLm6BI7WeqSqoq3iTkwQcxREZqvXAnlET11fsyqcMkulXpxod7PmTGrzNU6VUXUAn8KpJjMgkM8cPg4/ihbXt2NvHTZ2Bq1Ahz27YOb/9qvj72Nc9ve56GkQ2Z3W02Qb5BLjluaeTjZ8Ac5KOWkileRefnR/iwoWTt24fP4cNOOYZRb+SDdh9wV/W7mP7bdD7a+5FK4k6mEvhVJMdkXuypObzt5cuxxsYSMdI1ve/Pj37Oaztfo1V0K2Z0mYHZ6Pq15qVNcDm1qYnifYLvuQdjhQoEfLvSaYlVr9PzVuu3eLD2g8w/NJ/3fnlPJXEnMrg7gNJGSklybBa1m0c5vG17Rgbxs2bj36IF5ha3Obz9Ky05soT3fnmP9hXbM6HDBHz1vk4/picIjvTnxL6S1YguCWG3QkYCZCdDdgrkpEFuBuSmax+tOWDN1j7acsFu1R7SfkVDOtDpQWfQHgZfMPhpH43+4GMGnwDt4Rd06WE0gQfu6a6UjPDxIfyppzj/0kukb9xIYJcuTjmOTuh46baXMOqNLDq8CKvdysstXi6T1R3dTSXwK2Sl5ZGbZSXICWvAE5csxZaQQIQDKyMVZsGhBXy450M6VerEh+0/LHN1zUsiKNKf7Iw8sjPy8DM76LxICZmJkHoGUs5C2nlIj4G0C5AeC5nxkBEPmQm0z0mFbTfYrtDlJ2ij9u/LE6+0X0ruduuNx6r3Bf8wMIeBfzgEREJAFASW0z5aKkBQBQgsD+r3pkwJ6t2L05MnETd5CgGdOiF0zkmqQgiea/ocRp2RuQfnYpM2Xmv5mkriDqYS+BVSnDQD3ZaeQeKcOZjbtsW/cSOHtn2luQfnMnHvRLpW6cr77d7HqFNvwpe7OBM9NpNy1YowH8CaA0l/QcIJSDyh/Tvpb0j+G5JPg/XKpWkCzOFaUjSHQ3AVMIdzKiaFajc3utQj9g281FM2mrRHQU9ad4PzMKTUeut5WVqceZn5vfoMyE3TevrZKZCVDFlJ2sVGwUVF4glIiwFbzhXh67QkHlwFQqpASFUIvQlCq0PYTWByz8Y7SvEJg4GMnj0xzJ1H2g8/YOnuvOqLQghGNR6FQWdg1u+zyLPn8Vart9Df6O+0cl0qgV+hYAmZo++BJy1Zgi0lhYgRwx3a7pUKknePqj14p+07GHTqR3yl4PyfbUpMIQk8LwvijkLsYYj7Q/t33FEtUV8+jO0XpCW3iNpQsxsEVfx379UcCfr/nv+/t2yh2m0dHPufEiJ/CL2Yt0mk1Ib00y5oIwgFIwkpp7WLlFPb4LdlwGX3M/3DILy29v+PqA2RN0NUPe1iRSm1sps2JWLrNuKnTSOwWzen9cJBS+IjGo3AoDMw/dfpACqJO5B6d79CcmwWOp3AEubnsDZt6ekkzp2LuX07TLfe6rB2rzT/4HyVvG+AJdyEENrPmsxEOP8rnP8dzv8GFw5oPdKCRK0zQnhNKN8A6t8PYTUgrLrWE/UPde9/xJGE0HrUphAtEV/N5SMQCcch4RjE/QmHvtaSfwFzJETVhfK3auetfEMIqQZOTBRKEeh0hA8byrkxz5K2bh2WHj2cfsihDYYCMP3X6eiEjjdbvamG0x1AvcNfISUmE0uECZ3ecb9cSYsXa73v4c7rfS84tIAJeydwe9XbVfIuTF4WnPsV/dm9BPpVInn7V7Dv7UtfD6oM5epD3bsh6haIrKsl6qv0or2SwfdSb/tyUmr3+WMPQ8wh7eOFA7BrOtjztNf4BkF0Q6jQRHtUaq7de1fcwtK9O/EzZhBX0At3QTXIoQ2GIqVkxm8z0Akdr7d8XSXxElLvTFdIjs26OMTqCLa0NBLmzSegQwdM9Z1Tc3zR4UV8uOdDulbpyntt31PJu0BGPPyzC/7eBad/1nrY+QklWDeOFKpA17e0XmK5W8tWj9qVhIDAKO1RveOl5625EHdEO+/n9sPZvbBzyqUJdyFVodJt2qNKa+3CQM2Odwmh1xMxbBhnR48hde1agnr2dMlxhzYYik3amPX7LARCTWwrIfVOfxlpl6TEZlKxjuMm5yQuWoQ9JYVwJ/W+l/2xjA92f0CXyl14v9373p28M+Lhr+3a/dq/d2r3r0GbEBbdGFo+pSWLik0JWp3M+Z3nka3uc1k1PK9j8MkfQm8Ajftrz+VlabcrzvyiXVSd2Ay/L9e+5h8GVVpB1bZQrR1E1FEJ3YkCu3fHd8YM4qdNx9K9u0t64UIIhjccjpSS2QdmY9QZeem2l9TfYDF58bv9f2Wk5GDNsztsBrotLY3E+QsI6NgRU726Dmnzcl8f+5pxP4+jQ8UOfND+A++bbZ6XBX//qCWBk1sg5qD2vE8gVG4Bt/bRenbRjbRkcpngyFzycmxkpuZiDlLr413GaILKt2kPRmjD74kn80dKdsJfO+BI/q5W5kgtkVfvpD0s3le735mETkf4U09xdtQzpH6/hqBed7rmuPkT2/Lsecw/NB9fvS9jmo5RSbwYVAK/TMEGF46agZ60ZAn21FTCn3rKIe1dbvXJ1by+83VaRbfiww4fek/yTjgBf66D4+vhrx+1pU96Xy0hdHoVqrXXEvZ17ltfnIkem6kSuDsJoU0KDKsOjR7RniuY9X5qm3ZhdnCF9nzkLVCjM8EZUWBrrdaoO0Bgt2741qxB/MyZWHre4dQZ6ZcTQjC6yWhybDksOLwAX4MvIxqNcMmxyxKVwC+TnL+NaLADirjY0jNInDefgPbtHd77Xv/3el7e8TJNyzVlUsdJZbvCmi1P650dXQN/rtV6awDhtaDZE1rPrEor8Cnaz6ygUE9ybBbRNdV65lIlpAqE9IPG/bQeesxBOLEJjm+En2bS0J4Hf3yg3W+vfYe2hE/NXygWodMR9uST2ox0J68L/8+xheCF5i+Qa8tl1u+z8NP7MejWQS47flmgEvhlUuKy0Bt0BASXPCEmL/sMW0oK4cOGOiCyS3ac3cHYbWOpF16PqZ2mYjI4p2a7W+VmwPENcOQ7OPaDtkRJ76sNp7YYBjW7ahOgSiAw1BedXlzc+10ppYTQVgaUqw+tR0JOGgdXTqOe7zltJObwtyD0ULkl1LkDbu4FwZXdHbVHsXTvTvzUacTPmOn0deFX0gkdr7Z4lWxbNlP2T8FkMPHILY+47PieTiXwy6TEakvIhK5k92LsmZkkzJ2HuU0bTA0aOCg62Buzl2c2P0PN4JpM7zIdf6Nz9it3i5y0/Dfkb+DYBq2qmSkEavfQelnVO4FvgMMOp9PrsISbLlbeUzyEbyDxES2gQwew2+H8fvjjezj6Pax7SXtEN4Kbe0Pd/9OWASrXJPR6woc+ybmxz5O+aZPTaqQXRq/T87/W/yPHmsP7u9/HbDRzd827XRqDp3JpAhdCzAXuBGKllPVceewbkRLnmCVkScs/x5aY6NDe96GEQzy18SnKB5RnZteZWHwsDmvbbXIz4dg6OPiV1tO2ZkNAOWjUV3sDrtLaqWuwgyJNJMepHrjH0ukurSvv/Ko2P+LISji8Eja+qT3KN4R690DdeyC4krsjLrUsd9xB3NRpxE2fTkDnzi6fUGbQGXi/3fuM2DSCN3a9gdloplvVbi6NwRO5egHefMB1N1mKQNolKXFZJd7ExJ6dTcLcOfi3aIF/48YOie1E8gmeXP8kQT5BzOo6i1A/D77fZ7Nq9zK/fhI+rAlfDIR/foJG/eDRNTD6CPScADe1d3oBlaAIEylxWWq7w7IirDq0eQYGb4ZRB6Dr29oQ/PrXYFI9mNsD9szVqu8p/yIMBsKHDCHn8BHSt251Sww+eh8mdpjIreG38vz259lxdodb4vAkLk3gUsptQKn860lPzsGWZycoomQ98OQvVmCLi3dY7/tM2hkG/zAYg87A7G6zKWcu55B2XS7mEKx7GSbeAovv0YY96/4f9P8WxvwBPT/UJqO58P5bcKQ/1hwbmSm5Ljum4iLBlaH10zB4Czy9Hzq+om3e8t0z8GEtWNYX/litTZJUAG2nMmOFCsRPn+G2i1p/oz/TukyjRnANntn8DPti9rklDk8hXP2DEkJUBb671hC6EGIwMBggKiqqybJlyxx2/PT0dAIC/nsvNT1G8vdmSZWOgoCoYg4fWa2Ev/oatrAwkp4dU8JIIc2WxsQLE8mwZzAyaiTRPtElbrOkCjt/V2PISyMqZivlLmwiMP0EdqEnIawpMVEdSAxtil3vc/1GnCj9vOTvrZKqnQTmSNcNGRblHCr/VezzJyUB6SeJitlCVMw2fPKSyTVaiI1sz/nyXcgIqOrwWEuja50/0/btWJYsJenpp8m9pZCa+C6QZktj0oVJpNnSGFluJBV8Krgtliu54++3Y8eOe6WUTf/zBSmlSx9AVeDgjb6+SZMm0pE2b9581ecPbD0jpw7ZKFMTsordduLnn8vDtevItG3bi91GgbScNHn/yvtl00VN5f6Y/SVuz1EKO38X2WxSntgs5RePSflWhJSvW6Sc0UbKXTOkTI9zQYQ3Ljk2U04dslEe2nHWpce97jlUrskh58+aJ+XRtVIu7yflW+Ha7+knHaTcPUfKrJSSt1+KXev82XJy5J/t2su/HunnuoAKcTbtrOy0vJPsuLyjPJ162t3hXOSOv19gj7xKflRFaPOlxGaiNxZ/CZm0WkmY/Sl+detibtO6RLHk2HIYuXkkx5KO8VGHj2gY2bBE7blERjzsmAQfN4aFd2mFVpoMgCHb4cnt0OLJUrfN5KWlZGomutfRG6DW7fDAQhhzFLq/p+229t0zMKEOrByh1W/3MjofH8Iee5TM3bvJ3Ofe4evogGg+6foJObYchqwfQnxWvFvjKY1UAs+XEpelbTNZzCVkqWvWkvfPP4Q9OaREMzhtdhsvbHuBXy78wttt3qZtxbbFbsvppNQmoK14HD66GTa8ru2Dfc9sGPMn3DFe21KylCpYSpas1oJ7N/9QaDEUhv4IT2zSZq0fWAGzOsAn7WH/Yq1sr5cIvv9+9CEhxM+c6e5QqBFSg2mdpxGbGcuwDcNIz013d0iliksTuBDiM2AXUFsIcUYI8bgrj38tJVlCJu12EmbNwqdGdQI7dy52DFJK3vn5HTb8s4Hnmz3PnTe5pjZxkeVmwr6F8ElbmHs7HFsPTR+DYT/DY2vg1gfA6Lj91J2pYCa6oiAEVGwCd03VJlbe8aHWK//2Ke0C9YdXtf3Qyzidvz+hA/qTsW072YcPuzscGkY25KMOH3Es6RijNo8i16YmnRZw9Sz0h6SU5aWURillRSnlHFcevzAXl5AVcwZ6+ubN5Bw7RvjgwSWqYjTr91l8/ufnPFbvsdJZjSj1HNVOLtJmkq8coRXSuHMSjDkCPd6HyDrujrDI1FIy5ar8gqD5IBi2CwZ8p1UB3DUNpjTSZrD/9aM2AlVGhTz8MLqAAOI/meXuUABoW7Etb7V+i58v/MwrO17BLu3uDqlUUJXYgIyUXG0JWTHWgEspiZ/5CcZKlbDccUexY/j62NdM/XUqvW7qxajGo4rdjlOc2w87p8Lhb6hst0GdntqQY5XWHr/dY1CkCWuOjay0PPwt7p0Vr5RCQkC1ttoj5SzsmaOtJf/jO20P+ZbDtSH3Mraxit5iIaRvXxJmzSLnxAl8q1d3d0j0qt6LmMwYJu+bTKR/JM82e9bdIbmdugcOpMTl70JWjB545q5dZB84QNgTTyAMxbse2nZmG2/uepNW0a14s/WbpWNbPbsd/vwB5t+p3Qv8cx00H8LPt82EB5dA1TYen7wBgiK0izY1kU25rqAK0Pk1eOYw9JqsDa9/PRgmN4Afp0B2irsjdKjQAf0Rfn4kzJrt7lAuerze4zxU5yEWHF7AwkML3R2O26kEDhfvgRYngcfPno0hIoKgu/+vWMc+GH+QZ7c+S62QWnzU4SP3bwtqy4PflsGMVrD0fm33r27/g9GHoPs7ZJs8tJBMIQp+5uo+uHLDfPyhyUAY9hM8/IVWb339qzCxHqx/HdIuuDtChzCEhhJ8/32krF5N3rlz7g4H0HYwe77Z83St0pXxe8az5tQad4fkViqBAymxWej0goCQoi0hyzpwgMxdPxE6cCA6n6IPv55JO8NTG58i1C+U6V2mYzaai9yGw+Rmws+faPf4vh6i9a7vngUjf4NWI7R7gmVQYJgfQqgErhSDTge1usHA77SKbzW6wM4pMKk+rBp5aetbDxY2cCAACfPnuzWOy+l1et5t+y6NIxvz8o6X2XNhj7tDchuVwNGG0C3hJnT6op2OhNmfogsMJLjPA0U/Zk4KQzcMxWq3Mr3LdMJNblojnZOmrd+efCusGQtBFeHhz2HoTmjQp8zd27uS3qAjMMxPDaErJRPdCO6fB8P3QMO+8OtS+LgJfDUY4o66O7piM0ZHE9SzJ8lfrMCalOTucC7y1fsypdMUKgRUYOTmkZxM8fyLpeJQCRyKNQM959Qp0tavJ+Thh9EXsaxeri2XkZtHcjb9LFM6TeGmIDdseZidAlvHa72FDa9r+y0P/B4eW6sVuCgD97dvlFpKpjhMWHXoNUnbTKXFMDiyCqbdBp8P0PYD8EBhTzyOzMoiafESd4fyL0G+QczoMgODzsCwDcO8stCL1ydwKYu3hCxx7lyEjw+h/Yq23Msu7bzy4yvsjdnLuDbjaBLVpEjfX2LZqfmJ+1bY/D+o1EIrXtHva6hasgpynioowl8lcMWxAsvB7eNg1EFoO0bbgW9Gq/x0D99MAAAgAElEQVRE7v611UXhW7MmAZ06kbR4MfaMDHeH8y8VAysyrfM0ErISeHrT02RZvevv2OsTeFZaHnnZNixFSOB5MTEkf/MtwffegyG8aEPfU/dPZc2pNYxsPJIe1XoUNdziy0mH7RO0Hvfm/2k7fw3eCg8v04pXeDFLhImcTCvZ6WpnKsXBzGHaXuWjfod2z+Un8pbaNrpxf7o7uhsWNugJbCkpJK9Y4e5Q/qNeeD3eb/c+B+MP8sK2F7DZbe4OyWW8PoEXZwZ64oKFYLcT+thjRTrWN8e/YfaB2dxb814er+eiInR52bBrOkxpCBvfgkq3aRNuHvoMoj2gxroLqJnoitP5h0KnV7RE3vZZrXrh9Nvg66GQeMrd0V2Xf6NG+DdtSsK8+cjc0lcJrVPlTjzf/Hk2nd7EpH2T3B2Oy3h9Ak8t4hpwW2oqycuWYeneHZ+KFW/4OLsv7ObNXW/SonwLXm7xsvPXetusWrnTjxvDuhch8hZ4fAP0/VybcKNcVPCzT41XCVxxMv9QrUc+8jftHvmhr2BqU1g9BtJi3B3dNYUNHoT1wgVSVn/v7lCuqu/NfXmozkPMPzSfL/78wt3huITXJ/CUuCwQYAm7sQSetHw59sxMwp648R70Xyl/MWrzKCoHVmZChwnOXestJRxeqQ3TrRyhbS4yYBUMWAmVmjnvuB7MonrgiquZw7V75E//Co37w5552ijZpv+V2oIw5rZt8a1Zk8S5c0tt6eGxzcbSpkIbxv00jl3ndrk7HKdTCTw+i4BgX/TG658Ke24uSQsXYW7VCr+bb2yz++TsZJ7a+BR6oWdq56lYfCwlDblw//wEc7rC5/20z/sshic2aHWclUIZffT4B/mQonrgiqtZysOdE2H4bqjVHbaNh8kN4aeZYC1dQ9VCCEIff4ycY8fI2L7d3eFclUFnYHy78VQLqsaYLWM4kXzC3SE5ldcn8NQizEBPXbUKa1wcoY/f2L3vPFsez2x5hgsZF5jSaQqVAiuVJNTCxR/XNliYezuknIHeH8PQXXBzL69aDlYSQREmUlUPXHGXsOraOvLBW6BcPVj7vHaP/NA3pWrTlKA77sAQFUXCp6ViH6qrCvAJYFrnafjofXhq41MkZie6OySn8foEnhKXdUMz0KXdTsLcefjefDPmVq2u/3opGffzOPbE7OHN1m/SMNIJE8YyE+H7sdof+skt2iSZEXu1ITm92qemKILC1VpwpRSIbgT9V0LfFWDwgy8GaBfmZ/a6OzIAbensgAFk/vILWQcOujucQkUHRPNxp4+Jy4xj9JbR5NnK5goTr07gudlWstLybqgHnr51K7knThD22GM3NAFtyZElfHnsSwbVH+T4fb2tufkzyxvB7tlawn56v7ZMxceN5Vg9mCXCREZyDtZc71mCopRSQkDNrvDkDm00LfEUfNoJvhykjbC5WfAD96MLCCBhbunthQPUj6jPW63f0mpu/Dyu1N63LwmvTuAFs44t4ddP4Alz5mCILo+l++3Xfe2PZ39k/J7xdKrUieGNhpc4zn85tl4rCLHuRajQGJ78UbuHFhDp2ON4mUsz0bPdHImi5NPp8y/O92nFYA5/Cx83hS3vaXsXuIk+IICQB/uQtu4Hck+fdlscN6LnTT0ZVH8QXx77kiVHSlclOUfw6gR+o2vAs377jaw9ewkbOBBhvPYM8pMpJ3lu63PUDK7Ju23fRSccdIrjj8OSB2DJfYDU6pU/8hVE3eKY9r3cxZnoaiKbUtr4BmrbmI7YA7W7w5Z3YVpzOPiV2+6Ph/TrD3o9ifPmu+X4RTG80XA6V+7M+D3j2XF2h7vDcSiVwLl+Ak+YOw+dxULwvfdeu72cFEZsHIFRb2RKpyn4G/1LHmROurZF4fQW8PdO6Pq2NkHNy+qVO9vFHri6D66UVsGV4f75MHA1+AXDikdhQS+IPeLyUIxRkQT16kXyV1+Vqk1OrkYndLzT5h1qBtfkua3PcSql9BfOuVFencBT47PxNRvw9S+8V5175oy2aUmfB9CZC7+/bLVbGbttLOcyzjGp4ySiA6JLFpyU2hX2tObw4ySof782Qa3102Ao+talyrX5mY34+OlVD1wp/aq2gSFboecEuHAAZraBdS9r+xy4UOiAAcjsbJKXf+7S4xaHv9Gfjzt9jI/eh6c3PU1qrmvPlbN4dwKPyyToOve/ExcuBJ2OkEeuvWnJ5H2T2XluJy/f9jKNIktY6Sz+GCy8S7vC9g+Fx36Au2dAYFTJ2lUKJYTAopaSKZ5Cp4dmT8CIfdr2pbumwdRmcGCFy4bV/WrXwty6NYlLFmMvheVVr1Q+oDwT2k/gTNqZMlMz3asTeEp89jWXkNlSU0lZ8SVBPe/AGFV48lx1YhXzD83nwdoPcl+t+4ofUF4WbBqnTVI7tx96jNc2HKl8W/HbVG5YULhJlVNVPIs5DHpPgSc2ajugffk4LPo/SHBNAZPQRx/FFhdPaiktr3qlpuWa8uJtL7L97HY+3v+xu8MpMa9N4HabnfSE7GvOQE/+/HPsmZmEDhxY6GsOxR/ijZ1v0DSqKWObjy1+QCc2wfSWsO0DuOUuGL4HbhusXWkrLmEJN5Ean420l73lJkoZV7EJDNqkXfSf3ae9l2x5D6w5Tj2suXUrrbzqvHkes0zrgdoPcH+t+5lzcA5rTq1xdzgl4rUJPD0pB7tdFjqELvPySFy0GP8WLQotmxqfFc/IzSMJN4UXv8Z5epy2vnPR3Vqy7v8t3PupGi53A0uECZvVTkaKc9/0FMUpdHrton/4brj5Tm22+sw28NePTjukEILQgQPJ+fNPMnbudNpxHO3F5i/SOLIxr/34Gn8k/uHucIrNaxN4wWSlwobQU9euxRoTQ9ijA6/69Tx7HmO2jCElJ4XJnSYT6hdatACkhP1LYFozOPQ1tH9eW9N9U4eitaM4TMHFnBpGVzxaYDm4by70/RKs2TD/Dm1joyznzBa39LoTfUQ4ifMXOKV9ZzDqjUzoMIEg3yBGbR5Fcnayu0MqFq9N4AWTlSzhfv/5mpSSxHnz8aleHXPbtlf9/gl7JrAvdh+vt3qdOqF1inbwxFPafapvh0F4ba3iUseXwPjfWBTXsURo5z8lThVzUcqAml1g2E/Q6un8zsJt2k6FDqbz8SG0b18ytm8n59gxh7fvLOGmcCZ2mEhsZixjt431yElt3pvA47PR6QUBIf9Nmpm7d5N9+DChA/ojdP89RatOrGLJkSU8cvMjRSuTardpJVBntIIze7RlII+ugcgiXgAoThEQ6ofQCdUDV8oOHzN0e1u7Px4QCZ/3o+7B9xy+93hwnz4IPz9t1Y4HqR9Rn1davMKu87uYsn+Ku8MpMq9N4ClxWQSG+aHT/bcYSuLCheiDgwnq3fs/XzuScIQ3d71J06imjG46+sYPGPentinBuhehalt46mdtGchVLhAU99DrdQSG+qpNTZSyJ7ohDNoMnV8nLGGPVl/it2UOW3JmCAkh6K67SPl2JdZEz9r9656a93B/rfuZe3AuP/z1g7vDKRKvzR6p8VlXncCWe/o06Rs3EfxgH3R+/+6dJ2cnM2rzKIJ9g/mw/Yc3NmnNZoUdE7XJJAnH4Z7Z8PByCKroqP+K4kAWtZRMKav0Rmg7mt3NJkFEHfh6CHz2IKSec0jzof0eQebmkvx56S/scqUXmr/ArRG38sqPr3A86bi7w7lhXp3AA6+SwBMXLQKDgZCHHv7X8za7jRe2v0BcVhwTO0wkzBR2/YPEHYW53WDDG9ruQsN+hlsfUCVQSzGVwJWyLsu/Ijz6Pdz+LpzcCtNawK9LS9wb961RA3ObNiQtWYr0gMIul/PR+zCxw0T8Df48s+UZ0nPT3R3SDfHKBJ6TmUdOpvU/PXBbejopX36FpUd3jFH/3t1r5u8z+fHcj7zQ/AXqR9S/9gHsNtg5FWa21Sas3TsH+ixWS8M8gCXcj6y0PHKzre4ORVGcR6eHlsNg6I/ahkjfDIVlfSE9tkTNhg4YgDUujtS1ax0UqOtE+kcyvv14Tqed5rWdr3nEunavTOAFW0ZeOQM9ecUK7BkZhPYf8K/nt53ZxszfZtK7em/ur3X/tRtPPAXz74QfXoYanbVZoPXvU71uD1FQ2CctQc1EV7xAWHVtc5Ru4+D4Bm2m+qGvi92cuU1rfKpXJ3H+Ao9IgFdqVq4ZoxqPYv3f61l4uPRPyPPSBP7ffcClzUbSosWYmjTBVK/uxefPpJ3hxe0vUjukNq+0eAVRWCKWEvYugBmtIeYg/N8MeHCp6nV7GItaC654G50eWg2HJ7dDSBX4YqBWXCqr6GujhRCE9utH9uHDZO3d6/hYXWBA3QF0qdyFiXsnsufCHneHc01emcCvVsQlbdMm8s6eJXRA/4vP5dhyGL1lNFJKJnaYiMlQSNnV9Fj47CFY9TRUaAxDd0LDh1Wv2wMVjMoUjNIoiteIqA2Pb4AOL8HBL7Xlrie3FrmZoLt6ow8KInGB5xR2uZwQgrdbv02lwEo8t+054jLj3B1SobwygafFZ+Prb8DXZLj4XNLCRRijowns3Pnic+/98h5HEo8wrs04KlkqXb2xo2u0usMnNmmTQvqvhOBCXquUen5mI0Y/veqBK95Jb4AOz8MT68FogoW9Ye1LRaqprjOZCO7Th7SNm8g9c9aJwTpPgE8AEztMJCMvg+e2PYfVXjrnxHhlAk+Nz/rX8Hn2H3+QuXs3IX37IvTa5iGrTqxixZ8reKzeY3Ss3PG/jeRmwnejtWUYgeW1/XlbDlPruj2cEAJLmJqJrni5Ck1gyHatVsVP02B2J4g9csPfHvLQgyAESZ8tdWKQzlUjpAavtniVvTF7mbp/qrvDuSqvzDapCdn/msCWuHgxwmQi+L57ATiRfIK3f3qbJlFNGNFoxH8bOP87zOoAe+ZAy+EwaCNEXn3DE8XzWML9SFWT2BRv5+OvVYt8aDmkXdDe836ZfUPLzYzlyxPYtSvJX6zAnpnp/FidpFf1XtxX6z7mHJzD1tNFv53gbF6XwKVdkppwqQduTUoiddV3BPXW7ttk5mUyestoTAYTH7T7AIPOcNk3S/hpBnzaGbJToN83cPs4MPi66X+jOIMl3ERqXJZHzqJVFIer3R2G7dIqSH7/rDbfJyPhut8W2u8R7KmppKz6zgVBOs8LzV+gTmgdXtrxEufSHVP0xlG8LoFbs8BulRcTePIXK5A5OYQ+0hcpJW/99BanUk7xfrv3ifS/bC14Rjws7QNrX4DqnbWJatWvMrSueDxLuAlrnp3MVM8qRqEoThMQCX2/gO7vwYmNMLM1nNp2zW8xNW6M7y03k7R4kUdfDPvqffmo/UfYpZ1ntz6LVZae++Fel8BzM7SPlnA/pNVK0tKl+LdsgW/Nmnx57EtWn1zNsIbDaFG+xaVvOrVNWx52cgv0GA8PfQbmG6jEpngkNRNdUa5CCGgxFJ7YoG2SsqA3bPqfVi76qi8XhPbrT86x42T+/LOLg3WsSpZKvN36bQ7EH+CbpG/cHc5F3pfA8yvkWcJMpG3YiPXCBUL79eNo4lHe++U9WpZvyeBbB2svsttg87vaL6qfRbvXfdtgtTysjLOEqbXgilKo8g1g8FZo2Be2jddmqhdST91yRw/0oaEkLlrs4iAdr0uVLjxy8yNsTdvKxr83ujscAAzXf0nZkpchQUBgqB+nFy/CWLEiulbNeG5tXwJ9Anmn7TvohA5Sz8NXg+Cv7dDgYej5oXbV6SGsNjtZeTay8+zkWG3kWu3k2SR5NjtWu8Rmt2OXYLf/e2hLpxPoBOh1Ogw6gVGvw6gX+Bh0+Br0ZFklVpsdg77sXvsF5vfA0xJUAleUq/INgP+bBtXawXfPaCOUd38Ctbr962U6X1+C+zxAwsxPyD1zBp+Knr2J0+gmo9l2Yhuv7nyVOmF1qBBQwa3xeF0Cz82AgGBf8o4fJWvPXiLHjuWdPe/xV8pffNrtU8JN4dqa7i8HQV6mVlGt4cPXb9gJcqw2EtJzSUjPJT4jh6SMXBIzcknJyiM5M4+UrDxSs/NIzcojPcdKeraVjFwbmblW8mxOvOe0YQ1GvcDfx4DZR0+An4FAPyOBfgYsfkaCTEaC/bWPoWYfQsw+hJt9CQvwISzAB1+D3nmxOYDRR4/J4qOG0BXlehr00YpXfTEQlt4PrUdCp9e09eT5Qh58kITZn5K0ZClRz491X6wOYNQbeTT8USbETWDs1rHM7z4fo/4GdqV0Eq9L4HkZEBTmR9LSpQg/P3Y29mPlbysZ2mAozaOaaEPmW9/Xttt7YIFWncgJ0nOsnEnK5FxyFueSs7mQks35lGxiUrVHbFoOKVl5V/1enQCLyUiwyYjFpCXOyEA/zL4GzL56/H0MmIx6TD46TEY9vgY9RoPAR6/HoBcY9QK9TodOgE4ICm4ISMAuJXap9eCtdonVJsm12cizSnKsNg7+cYzoSlXJytMuFDJybGTkWEnLySMhPZe/4jNIztIuKuyFXEME+xuJCPAlyuJHlMWP8kF+lAvyo0KwiehgExVCTAT4uvdX0xLmR6rqgSvK9YXX1O6Lr30RfpwMp3fDfXPAEg2AMSqKwK5dSP7qKyKeHoHOVEhFSw8RbgznjVZv8OzWZ5myfwpjmo5xWyxel8Bz0yGwkp6Uxd8hbm/P24cm0qxcM4ZUvwcW36NNVGvwkLb+sQRD5lJK4tJyOBmfwV/xGfydmMk/CZn8k5jJ6aRMkjP/nZz1OkFkoC/lgvy4KcJMi5vCiAj0JSLQlzCzD2EB2scQfx8C/QzodO65D78l5y86dKh53dfZ7ZLU7DySMvNIzMglIT2HhIxc4tJyiEvLITYtm5jUHE6ciCc2LQfbFdk+2N9IpRB/Kof6UznMn6ph/lQJM3NTuJmIQN/Ca9I7iCXcxIWTKU49hqKUGUYT9JoElVvCd6O0nRjv/fTiSp3Qvn1JW7OW1NWrCb7vPjcHW3K3V72d3Rd2M//QfJqVa0a7iu3cEodXJXBbnh1rFhgvnERmZzOl2p+YDCbeq94H/awOkJUEvT+GRv1ueKKa1Wbn78RMjsWkcTw2neOx6ZyIy+BUfAbpOZdmZxr1gooh/lQK9efWikFUDPGnYkh+jzPYRESgL3o3JWVn0OkEwf4+BPv7UC382hdCVpuduPQcziVncTY5m7NJWZxJ0i52Dp1L4YfDF/51SyDA10C1cDPVI8zUiAygRmQgtaICqBJmdtg5tIT5cXxvLHabHV0Zvt+vKA7VoI82ye2LAbDobuj0MrQZg6lJE3xr1SJxyVKC7r3X6RfgrvBcs+f4NfZXXt7xMit6rSDK7PqNq7wqgaclavc0xa87iKsVwXbTaaaH9yJy6cMQVBEeXw/lby30+5MycjlyPpXD51M5cj6NPy6kciw2nVyr/eJrooP8qB4ZwH1NKnJThJmqYWaqhZuJDjaVqQTtSAa9jvJBJsoHmWhS5b9ft9kl55KzOBWfcfFxPDadn08l8s2vl2a/+hp01IwKoE45C3XKBXJLeQu3RFsI9vcpckyWcBPSLklPyvlX2V1FUa4jsg4M2gSrRmrLzM7sQdw9k5C+fbnw+utk7d+Pf+PG7o6yxHz1voxvP54+3/XhxR0vMrvrbPQ6187v8aoEXnBP03j2GPM7JNDfGE3b7dOgVg+4eyaYgi++NiE9h9/PpnDgTAoHzqZw6GwK51IuTWqKDPSldrlABraqSs3IAGpFBVIjMgCzm+/dlkV6naBSqDZ60a5WxL++lp5j5URsOkdj0jgWk8YfF9LY+mccK/aeufiaCsEm6kZbqF8hiPoVg7i1YjCh5msn9YKZ6FrZXZXAFaVIfMxwz2yo2BzWvQSftCeo96fEBgaStHhJmUjgANWCqvFi8xd5bedrzDk459ISZBdxebYRQnQHJgN64FMp5XuuOnbBrOI8fSLpNfSMOrYHOr1KdouRHDqfzv5/TrL/dDK/nU7mTFJWfrxQLdxMk6qhDIi2UDc6iJvLBxIWoMqnlgYBvgYaVAqmQaXgfz0fn57DkfOpHDqXysGzKRw+l8oPh2Mufr1SqIkGFYNpVDmERpWDqRtt+dfs+H+tBa8d4pr/jKKUJUJodTOiG8LnA9At7U1wuztIXPcDkbGxGCMjr9+GB/i/Gv/HrnO7mP7rdJqXa07DyIYuO7ZLE7gQQg9MA7oCZ4DdQoiVUsrDrjh+4okLCLuVH29O5q3ETFbUmsgXB2tycN0Gcm3aMHiFYBMNKwXTr0UVbq0YTL0KFgL93LdMQCme8ABf2taMoG3NSz321Ow8Dp1N5bczyfx+Jpl9fyfx3e/nAfDR66hfMYimVUJoUiWEJpWCEQLS1KYmilIylZpruzV+MZCQpG9ItEaRvHw5ESOuslGUBxJC8GrLV/k9/nfGbhvLit4rsPhYXHJsV/fAmwPHpZQnAYQQy4C7AOcncCk5tmMHvtZyRFUTDEx+jdjUKOpXFAxsXZXGlUNoXDmYSIvf9dtSPJLFz0jL6mG0rH6pDG5Majb7/0li79/aY96Pf/HJtpMAPKU38cvBWG7rfZO7QlaUsiEgEvp/i8/61zDvW0bygk8I73cfIri8uyNziECfQMa3G0//Nf15Y+cbTGg/wSUT9YQri8wLIe4Dukspn8j/vB9wm5Ry+BWvGwwMBoiKimqybNmyEh9bSsnJBdtB+nC8TV1uCvenikWHj15NLCuK9PR0AgIC3B2G0+TaJH+l2jmaaMPviMDgCx27OnYEpqyfQ2dT569k3H3+IrcvxXftdmKefpHcqMpui6O4rnX+NqRs4M/sP3ki4gl8dEWfPFuYjh077pVSNr3yeVcn8PuB269I4M2llIWOpTRt2lTu2bPHMQHYbWzavJlOnbs4pj0vtGXLFjp06ODuMDyaOoclo85fybj7/EkpwWZFGDzz1uS1zp9dardidcKxS0+FEFdN4K4eQj8DVLrs84qA6zZY1enR6dUscUVRFHcRQoCHJu/rcXTivu7xXHo02A3UFEJUE0L4AA8CK10cg6IoiqJ4PJd2R6WUViHEcGAd2jKyuVLKQ66MQVEURVHKApePJ0spvwe+d/VxFUVRFKUsUUWeFUVRFMUDqQSuKIqiKB5IJXBFURRF8UAqgSuKoiiKB1IJXFEURVE8kErgiqIoiuKBVAJXFEVRFA+kEriiKIqieCCXbmZSHEKIOOBvBzYZDsQ7sD1vo85fyalzWDLq/JWMOn8l447zV0VKGXHlk6U+gTuaEGLP1XZ1UW6MOn8lp85hyajzVzLq/JVMaTp/aghdURRFUTyQSuCKoiiK4oG8MYHPcncAHk6dv5JT57Bk1PkrGXX+SqbUnD+vuweuKIqiKGWBN/bAFUVRFMXjqQSuKIqiKB7IaxK4EKK7EOKoEOK4EOIFd8fjaYQQlYQQm4UQR4QQh4QQI90dkycSQuiFEPuFEN+5OxZPI4QIFkKsEEL8kf972NLdMXkaIcQz+X+/B4UQnwkh/NwdU2kmhJgrhIgVQhy87LlQIcR6IcSx/I8h7orPKxK4EEIPTAN6ALcADwkhbnFvVB7HCoyRUt4MtACeUuewWEYCR9wdhIeaDKyVUtYBGqDOY5EIISoATwNNpZT1AD3woHujKvXmA92veO4FYKOUsiawMf9zt/CKBA40B45LKU9KKXOBZcBdbo7Jo0gpz0sp9+X/Ow3tzbOCe6PyLEKIikBP4FN3x+JphBAWoB0wB0BKmSulTHZvVB7JAJiEEAbAHzjn5nhKNSnlNiDxiqfvAhbk/3sB8H8uDeoy3pLAKwCnL/v8DCr5FJsQoirQCPjZvZF4nEnAWMDu7kA80E1AHDAv/xbEp0IIs7uD8iRSyrPAh8A/wHkgRUr5g3uj8khRUsrzoHVsgEh3BeItCVxc5Tm1fq4YhBABwJfAKCllqrvj8RRCiDuBWCnlXnfH4qEMQGNghpSyEZCBG4cuPVH+vdq7gGpANGAWQjzi3qiUkvCWBH4GqHTZ5xVRQ0dFJoQwoiXvJVLKr9wdj4dpDfQWQvyFdgunkxBisXtD8ihngDNSyoJRnxVoCV25cV2AU1LKOCllHvAV0MrNMXmiGCFEeYD8j7HuCsRbEvhuoKYQopoQwgdt4sZKN8fkUYQQAu3+4xEp5UfujsfTSClflFJWlFJWRfv92ySlVL2fGySlvACcFkLUzn+qM3DYjSF5on+AFkII//y/586oiYDFsRIYkP/vAcC37grE4K4Du5KU0iqEGA6sQ5t5OVdKecjNYXma1kA/4IAQ4tf8516SUn7vxpgU7zICWJJ/EX4SeNTN8XgUKeXPQogVwD60VSX7KUVlQUsjIcRnQAcgXAhxBngdeA/4XAjxONpF0f1ui0+VUlUURVEUz+MtQ+iKoiiKUqaoBK4oiqIoHkglcEVRFEXxQCqBK4qiKIoHUglcURRFUTyQSuCKUkYIIT4XQiQKIcpd8bxeCLEnf/ck0zW+/w0hhMx//OXAuF65rN0zjmpXUbydSuCKUnYMRysRPP2K559Fq1r2hJQy6wbaaQnc7cC45uW3qWoGKIoDqQSuKGWElDIWeAa4WwhxP4AQohbwBvCJlHLrDbbzk5RyvwPjOiul/AltMxJFURxEJXBFKUOklAuBtcBUIUQEWvnbOOD54rYphKiaP/z9pBDiXSHEBSFEmhBicX5ZzhpCiHVCiHQhxHEhxIDrt6ooSkmpBK4oZc8QtL2efwLaAE86aOe4F9F2sRoAvAb0AWYCXwOr0Ybdf0fb8rOuA46nKMo1eEUtdEXxJlLKf4QQU9G22/zKgfXqT0gpC3rX64QQbdHq4/eTUi4GEELsAXoD9wFqvwFFcSLVA1eUMkYIYUFLrBJoJoQIdFDTa674/I/8j+sKnpBSJqFtr1gJRVGcSiVwRSl7xgMhQE8gEnjXQe0mXfF57jWe93PQMRVFKYRK4IpShggh2gODgFeklGuA/wHDhBCt3BuZoq52SwUAAADxSURBVCiOphK4opQR+UVaPgV2A5Pzn34fOAh8mr+PtqIoZYRK4IpSdrwFVEEr2GIHkFLmAY8DtYGX3RiboigOphK4opQBQoimaEVc3pNSHrj8a1LKgh75C2p5l6KUHUJK6e4YFEUpBYQQbwCvA0ZASiltDmpXAHq0ojKdpZQVHdGuong71QNXFOVKecAJB7b3cn6b/R3YpqJ4PdUDVxQFACFENFqlNYCcK4fiS9BueaBC/qe5UsrfHdGuong7lcAVRVEUxQOpIXRFURRF8UAqgSuKoiiKB1IJXFEURVE8kErgiqIoiuKBVAJXFEVRFA/0/x6etASWbs4VAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# We are going to plot the trajectories, of a series of launchs that start\n",
    "# with different angles\n",
    "thetas = [10, 30, 45, 60, 80]\n",
    "\n",
    "# We will keep the initial speed fixed\n",
    "v_0 = 10\n",
    "\n",
    "# Compute, and plot, each trajectory\n",
    "font = 16   # set font size\n",
    "\n",
    "fig = plt.figure(figsize=(7, 5))\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.set_xlabel(r\"X [m]\", fontsize=font), ax.set_ylabel(r\"Y [m]\", fontsize=font)\n",
    "\n",
    "for theta in thetas:\n",
    "    v_t0 = initial_velocity(v_0, theta)\n",
    "\n",
    "    ts = np.linspace(0, flight_time(v_0, theta), 100)\n",
    "\n",
    "    xs, ys, zs = parabolic_positions(ts, (0, 0, 0), v_t0)\n",
    "\n",
    "    ax.plot(xs, ys, label=r\"{} deg\".format(theta))\n",
    "\n",
    "ax.legend(), ax.grid()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:37:32.927034Z",
     "start_time": "2021-01-31T23:37:32.716967Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ea64c59905f8406cbfa68e0bbe4d21ff",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=41, description='theta', max=89, min=1, step=5), Output(layout=Layout(he…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# As we did in the week 1 notebook, we can actually animate this\n",
    "# plot to, interactively, see the influence of the angle on the reach\n",
    "\n",
    "from ipywidgets import interactive       # make interactive plots\n",
    "\n",
    "# These properties are common to all the plots\n",
    "r_0 = (0, 0, 0)\n",
    "v_0 = 10\n",
    "\n",
    "# Define the values for the parameters we want to visualize\n",
    "theta_min, theta_max = 1, 89\n",
    "theta_step = 5\n",
    "\n",
    "\n",
    "def f(theta):\n",
    "    # initialize velocity\n",
    "    v_t0 = initial_velocity(v_0, theta)\n",
    "    # compute the flight time\n",
    "    t_f = flight_time(v_0, theta)\n",
    "    ts = np.linspace(0, t_f, 100)\n",
    "    # compute the ball's positions\n",
    "    xs, ys, zs = parabolic_positions(ts, r_0, v_t0)\n",
    "    # compute the maximum distances for scaling\n",
    "    x_max = v_0**2 / g\n",
    "    y_max = v_0**2 / g / 2\n",
    "    # plot\n",
    "    plt.rcParams.update({'font.size': 14})\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(1, 1, 1)\n",
    "    ax.plot(xs, ys)\n",
    "    ax.set_xlim((0, x_max)), ax.set_ylim((0, y_max))\n",
    "    ax.set_xlabel(r\"X [m]\"), ax.set_ylabel(r\"Y [m]\")\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "# Output interactive plot\n",
    "interactive_plot = interactive(f, theta=(theta_min, theta_max, theta_step))\n",
    "output = interactive_plot.children[-1]\n",
    "output.layout.height = '350px'\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## b. Numeric approach\n",
    "\n",
    "In general, the differential equations of motion do not have an analytical solution, but we can compute one numerically. The method we are going to use here is not how it is usually done (there are far better methods), but it is the simplest to illustrate the principles behind the numerical solution of differential equations.\n",
    "\n",
    "In general, we want to solve equations of the form:\n",
    "\n",
    "$$\\frac{dy}{dt} = f(y, t)$$\n",
    "\n",
    "We can start with a known initial condition:\n",
    "\n",
    "$$y(t_0) = y_0$$\n",
    "\n",
    "Then, the solution for another time, close to the initial time, can be approximated as:\n",
    "\n",
    "$$y(t_1) \\approx y(t_0) + \\Delta y |_0 = y(t_0) + \\frac{\\Delta y}{\\Delta t}|_0  \\Delta t= y_0 + f(y_0, t_0) \\Delta t$$\n",
    "\n",
    "And iteratively, we can compute the solution at any point in time using the solution at a previous point in time, and evaluating the function $f$ at that point in time:\n",
    "\n",
    "$$y(t_{i+1}) \\approx  y_i + f(y_i, t_i) \\Delta t$$\n",
    "\n",
    "Let's implement such a function, and use it to solve our problem:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can decompose our problem in the following steps:\n",
    "\n",
    "\n",
    "  1. Using the velocity, compute the new position at the next point\n",
    "  $$\\mathbf{x}_i = \\mathbf{x}_{i-1} + \\mathbf{v}_{i-1} * \\Delta t$$\n",
    "\n",
    "  2. Using the acceleration, compute the new velocity\n",
    "  $$\\mathbf{v}_i = \\mathbf{v}_{i-1} + \\mathbf{a}_{i-1} * \\Delta t$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:37:47.369615Z",
     "start_time": "2021-01-31T23:37:47.322234Z"
    }
   },
   "outputs": [],
   "source": [
    "def parabolic_positions_numeric(t, r0, v0):\n",
    "\n",
    "    import numpy as np\n",
    "\n",
    "    # as before t is an array with times, let's initialize the outputs\n",
    "    xs = np.zeros(len(t))\n",
    "    ys = np.zeros(len(t))\n",
    "    zs = np.zeros(len(t))\n",
    "\n",
    "    # we start at the initial condition\n",
    "    xs[0] = r0[0]\n",
    "    ys[0] = r0[0]\n",
    "    zs[0] = zs[0]\n",
    "\n",
    "    # with the initial velocity\n",
    "    vx = v0[0]\n",
    "    vy = v0[1]\n",
    "    vz = v0[2]\n",
    "\n",
    "    # let's iterate to compute the positions at all subsequent points in time\n",
    "    for i in range(1, len(t)):\n",
    "        # compute the time step\n",
    "        dt = t[i] - t[i-1]\n",
    "        # update positions with velocity\n",
    "        xs[i] = xs[i-1] + vx * dt\n",
    "        ys[i] = ys[i-1] + vy * dt\n",
    "        zs[i] = zs[i-1] + vz * dt\n",
    "        # update velocity with acceleration\n",
    "        vx = vx\n",
    "        vy = vy - g * dt\n",
    "        vz = vz\n",
    "\n",
    "    return xs, ys, zs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:37:49.806862Z",
     "start_time": "2021-01-31T23:37:49.605327Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "24a0a7b2447d41e5a983ac29e8c6c517",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=50, description='time_steps', min=5, step=5), Output(layout=Layout(heigh…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# The precision of the approximation will depend on the time step used\n",
    "# Let's see that with an interactie plot\n",
    "\n",
    "# common parameters that will not change\n",
    "r_t0 = (0, 0, 0)\n",
    "v_t0 = initial_velocity(10, 45)\n",
    "t_f = flight_time(10, 45)\n",
    "\n",
    "# parameters that we want to change\n",
    "min_time_steps = 5\n",
    "max_time_steps = 100\n",
    "change_time_steps = 5\n",
    "\n",
    "\n",
    "def f(time_steps):\n",
    "    # compute the times that correspond to the number of time steps desired\n",
    "    ts_a = np.linspace(0, t_f, 100)\n",
    "    ts_n = np.linspace(0, t_f, time_steps)\n",
    "    # compute the positions using the exact, analytic formula\n",
    "    xs_a, ys_a, zs_a = parabolic_positions(ts_a, r_t0, v_t0)\n",
    "    # compute the positions using the approximate, numeric method\n",
    "    xs_n, ys_n, zs_n = parabolic_positions_numeric(ts_n, r_t0, v_t0)\n",
    "\n",
    "    # plot\n",
    "    plt.rcParams.update({'font.size': 14})\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(1, 1, 1)\n",
    "    ax.plot(xs_a, ys_a, 'k-', label='exact')\n",
    "    ax.plot(xs_n, ys_n, 'b-', label='approximate')\n",
    "    ax.legend()\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "# Output interactive plot\n",
    "interactive_plot = interactive(f, time_steps=(min_time_steps, max_time_steps,\n",
    "                                              change_time_steps))\n",
    "output = interactive_plot.children[-1]\n",
    "output.layout.height = '350px'\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see how the numeric approximation tends to overshoot the exact solution, and needs many time steps (around 500) to fit the known solution.\n",
    "\n",
    "This is because we are using a poor numeric approximation. For example, to propagate the position in y we are using the vertical position in the previous position, which is the maximum vertical velocity as the ball moves in that time interval. \n",
    "\n",
    "A simple improvement would be to update the velocity first, to the average between the velocity at the beginning and the end of the interval:\n",
    "\n",
    "  1. Using half the acceleration, compute the new velocity\n",
    "  $$\\mathbf{v}_i = \\mathbf{v}_{i-1} + 0.5 * \\mathbf{a}_{i-1} * \\Delta t$$\n",
    "\n",
    "  2. Using the velocity, compute the new position at the next point\n",
    "  $$\\mathbf{x}_i = \\mathbf{x}_{i-1} + \\mathbf{v}_{i-1} * \\Delta t$$\n",
    "\n",
    "  3. Using half the acceleration, update the new velocity\n",
    "  $$\\mathbf{v}_i = \\mathbf{v}_{i-1} + 0.5 * \\mathbf{a}_{i-1} * \\Delta t$$\n",
    "\n",
    "This scheme is called \"leap frog\" and it improves numerical convergence dramatically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:38:02.039835Z",
     "start_time": "2021-01-31T23:38:02.028367Z"
    }
   },
   "outputs": [],
   "source": [
    "def parabolic_positions_numeric(t, r0, v0):\n",
    "\n",
    "    import numpy as np\n",
    "\n",
    "    # as before t is an array with times, let's initialize the outputs\n",
    "    xs = np.zeros(len(t))\n",
    "    ys = np.zeros(len(t))\n",
    "    zs = np.zeros(len(t))\n",
    "\n",
    "    # we start at the initial condition\n",
    "    xs[0] = r0[0]\n",
    "    ys[0] = r0[0]\n",
    "    zs[0] = zs[0]\n",
    "\n",
    "    # with the initial velocity\n",
    "    vx = v0[0]\n",
    "    vy = v0[1]\n",
    "    vz = v0[2]\n",
    "\n",
    "    # let's iterate to compute the positions at all subsequent points in time\n",
    "    for i in range(1, len(t)):\n",
    "        # compute the time step\n",
    "        dt = t[i] - t[i-1]\n",
    "        # update velocity with half acceleration\n",
    "        vx = vx\n",
    "        vy = vy - g * dt / 2\n",
    "        vz = vz\n",
    "        # update positions with velocity\n",
    "        xs[i] = xs[i-1] + vx * dt\n",
    "        ys[i] = ys[i-1] + vy * dt\n",
    "        zs[i] = zs[i-1] + vz * dt\n",
    "        # update velocity again with half acceleration\n",
    "\n",
    "        # START YOUR CODE HERE (3 lines of code)>>>>>>>>>>>>>>>>>>\n",
    "        vx = vx\n",
    "        vy = vy - g * dt / 2\n",
    "        vz = vz\n",
    "        # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> END YOUR CODE HERE\n",
    "\n",
    "    return xs, ys, zs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:38:03.772649Z",
     "start_time": "2021-01-31T23:38:03.765319Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0., 1.]), array([0.    , 0.8038]), array([0., 0.]))"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's test the output for the flight_time helper function\n",
    "parabolic_positions_numeric([0, 0.2], (0, 0, 0), (5, 5, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Expected Output**: \n",
    "<table style = \"width:40%\">\n",
    "    <tr>\n",
    "    <td>** parabolic_positions_numeric([0, 0.2], (0,0,0), (5,5,0)) **</td> \n",
    "        <td>(array([0., 1.]), array([0.    , 0.8038]), array([0., 0.])) </td> \n",
    "    </tr>\n",
    "\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:38:05.754793Z",
     "start_time": "2021-01-31T23:38:05.573710Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "615e4deb5f9a44f2b1b693fbd35161ca",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=50, description='time_steps', min=5, step=5), Output(layout=Layout(heigh…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's see how the numeric precision with this improved method is so much\n",
    "# better\n",
    "\n",
    "# common parameters that will not change\n",
    "r_t0 = (0, 0, 0)\n",
    "v_t0 = initial_velocity(10, 45)\n",
    "t_f = flight_time(10, 45)\n",
    "\n",
    "# parameters that we want to change\n",
    "min_time_steps = 5\n",
    "max_time_steps = 100\n",
    "change_time_steps = 5\n",
    "\n",
    "\n",
    "def f(time_steps):\n",
    "    # compute the times that correspond to the number of time steps desired\n",
    "    ts_a = np.linspace(0, t_f, 100)\n",
    "    ts_n = np.linspace(0, t_f, time_steps)\n",
    "    # compute the positions using the exact, analytic formula\n",
    "    xs_a, ys_a, zs_a = parabolic_positions(ts_a, r_t0, v_t0)\n",
    "    # compute the positions using the approximate, numeric method\n",
    "    xs_n, ys_n, zs_n = parabolic_positions_numeric(ts_n, r_t0, v_t0)\n",
    "\n",
    "    # plot\n",
    "    plt.rcParams.update({'font.size': 14})\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(1, 1, 1)\n",
    "    ax.plot(xs_a, ys_a, 'k-', label='exact')\n",
    "    ax.plot(xs_n, ys_n, 'b-', label='approximate')\n",
    "    ax.legend()\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "# Output interactive plot\n",
    "interactive_plot = interactive(f, time_steps=(min_time_steps, max_time_steps,\n",
    "                                              change_time_steps))\n",
    "output = interactive_plot.children[-1]\n",
    "output.layout.height = '350px'\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See how the performance improves tremendously, and now only ~20 time steps suffice to get an answer that is essentially the same as the exact solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding air drag to our problem\n",
    "\n",
    "+ The big advantage of using a numerical approach, is that now we can solve problems that are more complex, and for which we do not have a good analytic solution.\n",
    "\n",
    "+ For example, let's remove the simplifying assumption about air drag. Let's now suppose that the ball experiences a drag that is proportional to its speed, and oposses it:\n",
    "\n",
    "$$\\mathbf{F_D} = - C_D \\mathbf{v}$$\n",
    "\n",
    "+ Now, we can modify the numeric function just changing the acceleration term:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:38:18.158170Z",
     "start_time": "2021-01-31T23:38:18.147589Z"
    }
   },
   "outputs": [],
   "source": [
    "def parabolic_positions_numeric_drag(t, r0, v0, c_d=0):\n",
    "\n",
    "    import numpy as np\n",
    "\n",
    "    g = 9.81\n",
    "\n",
    "    # as before t is an array with times, let's initialize the outputs\n",
    "    xs = np.zeros(len(t))\n",
    "    ys = np.zeros(len(t))\n",
    "    zs = np.zeros(len(t))\n",
    "\n",
    "    # we start at the initial condition\n",
    "    xs[0] = r0[0]\n",
    "    ys[0] = r0[0]\n",
    "    zs[0] = zs[0]\n",
    "\n",
    "    # with the initial velocity\n",
    "    vx = v0[0]\n",
    "    vy = v0[1]\n",
    "    vz = v0[2]\n",
    "\n",
    "    # let's iterate to compute the positions at all subsequent points in time\n",
    "    for i in range(1, len(t)):\n",
    "        # compute the time step\n",
    "        dt = t[i] - t[i-1]\n",
    "        # update velocity with half acceleration\n",
    "        vx = vx - c_d * vx * dt / 2\n",
    "        vy = vy - c_d * vy * dt / 2 - g * dt / 2\n",
    "        vz = vz - c_d * vz * dt / 2\n",
    "        # update positions with velocity\n",
    "        xs[i] = xs[i-1] + vx * dt\n",
    "        ys[i] = ys[i-1] + vy * dt\n",
    "        zs[i] = zs[i-1] + vz * dt\n",
    "        # update velocity again with half acceleration\n",
    "        vx = vx - c_d * vx * dt / 2\n",
    "        vy = vy - c_d * vy * dt / 2 - g * dt / 2\n",
    "        vz = vz - c_d * vz * dt / 2\n",
    "\n",
    "    return xs, ys, zs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:38:19.086185Z",
     "start_time": "2021-01-31T23:38:18.937170Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEQCAYAAAD2/KAsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3xUVfrH8c+ZhEkCIbQAoSb03gMYalRY1BVdK6IioCgKNhAQQcACsiqgwEqz74LC/tC1wSIs0lsIAWmhSAktEHpJISTz/P6YGANSAiQ5k8zzfr3mBZl779wvl8k8c+499xwjIiillFJ5zWE7gFJKKe+kBUgppZQVWoCUUkpZoQVIKaWUFVqAlFJKWaEFSCmllBW+tgNcKjg4WMLCwmzHUEoplQPWrVt3TERKX26ZxxWgsLAwoqOjbcdQSimVA4wxcVdapqfglFJKWaEFSCmllBVagJRSSlmhBUgppZQVWoCUUkpZoQVIKaWUFVqAlFJKWaEFSCmllBUedyOqUvnViRPJ/PLLflauPM62bec5dMhw7JiTM2eKkJoaQFpaAC5XEUSy/toJDkcSvr5J+PqmUKRIIiVLJlO2rIuwMB/Cw4O47bby1KkTjMNhrP3blMoN5lozohpjXgPuB2oB54HVwGsisvkq24QBey6z6E4RmXe1/YWHh4uOhKA83bFjSUyfvoP//e8UGzc6iY+vQFpaJS4+qXAGp/MohQufISDgPAEBaRQpko7T+cfvXHq6ISnJQXKyD8nJhUhMDOT8+ZK4XKUvei1jTlKs2F5q1DhNmzZ+dOkSSsuW5fPs36vUjTLGrBOR8Msuy0YB+hmYCawFDPAWEAHUFZETV9gmDHcBugP4NcuiEyKSerX9aQFSnig1NZ3PP9/KzJnHWb++OKdP1wH8APDxOUBIyH5q1DhPgwZOWrcuRdu25SlfvugN7y8lJY2oqHiWLj1CTEwi27YZ9u0rRWJi9cz9+vrGUb16HB07+tKvX12qVCmeA/9SpXLWTRWgy7xYIHAa+JuI/HiFdcJwF6DmInJd1UQLkPIUZ86c5/33f2XWrPP89ltt3OMpuihceBv16iVw111F6NKlCnXqBOdppu++28V33x1l1Sp/Dh+uDRQD0ihRYiO33nqWQYNqaOtIeYycLkDlgENAWxFZfoV1wnAXoP2AP7AT+EBEZl/r9bUAKZvS0lxMnryJSZNOs317Q0SKA6cJDd3MPfcYj2tppKSkMX36Nj777BgxMRU5f7464KJ48V+5995zvPlmQ0JDi9mOqbxYThegfwM1gHARSb/COsFAd2AFkAbcAwwFuovI9Mus/wzwDEDlypWbxcVdcfBUpXLF9u3H6d9/E/PnVyEtLRRIpEqV9Tz5pD8vv9yQwECn7YjZsmDBXv7+970sWxbKhQtVgGSqVVvH0KEl6N69rnZkUHkuxwqQMWYc8AjQRkR2X2eISRnbNbzaetoCUnnp3//ezmuvHWX37nDAn6CgX3n00XOMGNGIkJBA2/FumMsl/Otfsbz77jFiY5sCgQQExPL448cZN65FvimoKv+7WgHK9n1AxpgPgK7AbddbfDKswd1yUsoql0sYN249wcHRdOlSi927G1O//hq+/XYnp083YvLk1vm6+AA4HIbu3euydWs79u930aXLElwuXz7+uA3Fix/jb39bzKFDZ23HVF4uWwXIGDMeeBR38dl2g/tqDMTf4LZK5YixY9dTrNhmXnmlCSdOhNKx42J2705j06b23Hdfwfx+VLFiEDNnticpqTpvvLGWwMAEvv8+kooV07jjjsUkJCTajqi81DULkDHmI6An7tbPSWNMSMYjMMs6o40xC7P83N0Y86gxpo4xppYxZgDQF5iYC/8Gpa5p8uSNlCixngEDmpCcXIpHHlnKiRNFmT8/0qM6FeQmh8MwYkRzTp1qzCefbKZ06Z38/HMkISFJ/O1vizl1KsV2ROVlstMC6gMUBRbibsH8/hiQZZ1yQLVLtnsdiMZ9/9AjwJMi8sHNBlbqeixevI+KFVfRp09DzpypwP33L+HYsZJ8/XU7ihf3tx3Pmqeeqs+RIy2YOnUTxYvv5/vvIyld+ij9+q3E5bq+jklK3ajr7gWX27QTgsoJBw6c4b77YoiOjgDSuO22tcya1YLg4MK2o3mk99+PYfjwIqSk1KJo0Y1MnOhL9+51bcdSBUCOdEJQKj9wuYR+/VYSGppEdHQ7atSIYt26syxcGKnF5yoGDmzK6dPV6d59GYmJIfToUZv69ZcSF3fadjRVgGkBUgXGL7/EUabMOj78sBV+fif54ott7NjRlqZNQ2xHyxecTh+++KItcXH+NGmyjC1bWlO16nn69l2hp+VUrtACpPK9tDQX99+/hNtvL83x4zV54IElnDpVU08h3aCKFYOIiWnPV1/tJCAggUmTWlOu3FpiYg7bjqYKGC1AKl9buDCOUqU28Z//tKd06S2sXZvI7NntcTp9bEfL97p2rc2JE3W4774lJCTUJzzcn6efXq6tIZVjtACpfMnlEh5/fCkdOgRz5kwYvXot5/DhcMLDy9mOVqA4nT58+2175s9PoGjROD75pA0VKkSxfftx29FUAaAFSOU7O3eeoGLFNcyY0Y6SJbexdm0SH3/cRsc5y0UdO4Zx/HgD7rlnMYcPN6Zu3Qu8/36M7Vgqn9MCpPKV8eM3ULv2eeLjm3L33Ys5cqSJtnryiK+vg++/j2TmzL34+iYyaFBjWrZcTFLSBdvRVD6lBUjlCy6XcOedi3n55Qb4+KQwffoufvwxEl9ffQvntS5danHwYDlq115OVFQk5cpt1Q4K6obob6/yeHFxp6lQIYp58yKpVGkNe/cG89hjdWzH8mrBwYWJjW1Hnz4rOHOmOs2bO/jwww22Y6l8RguQ8mjff/8b1auf4vDhptx33xL27o24qamuVc766KPWfPddPL6+5+jXrwH33rtYe8mpbNMCpDzWsGFR/O1vZXG5/JkyZRvfftteOxp4oHvvrc6ePaWpUCGKH36IpGbNFTqwqcoWLUDK47hcQqdOixk5MpzChfezerWL3r0b2I6lrqJ8+aLs3duSyMjF7NrVhkqVdrJhwxHbsZSH0wKkPMq5c6nUqLGC+fMjCQ1dxf79VWjeXHu55Qe+vg4WLYpkwIBVnDtXlfDwdL75ZoftWMqDaQFSHiMu7jSVKm1h9+423HbbYnbvbkXJkgG2Y6nr9P77EcyceQARw4MPluW99/R+IXV5WoCUR1i16iC1aiVw6lQ9evdezsKFkXq9Jx/r0qUWK1cK/v5HePXVBvTqtdx2JOWBtAAp6/7zn520aePg/PkyjBmzhSlT2tiOpHJAy5bl2b69LCVLbuTTT9vwl79oDzl1MS1AyqpPPtnMAw8EA4Zvv03glVea2I6kclDlysWIi2tAWNgKFiyIpEWLJVqEVCYtQMqa0aPX8fTTVfD1Pc2iRRe4774atiOpXBAY6GTnzggaNlzCunWR1Kq1nJSUNNuxlAfQAqSsGDRoNUOGNMDf/wAxMUVo166S7UgqF/n6Oli/vh2RkYv57be2VK0azblzqbZjKcu0AKk899JLK3n//XACA7ezbVsI9euXth1J5QGHw7BoUST337+E+PhbqFZtg96w6uW0AKk89dxzK5gwoSXFim1m+/YwQkOL2Y6k8tg337Tn0UeXkpDQgurVN3PiRLLtSMoSLUAqz/TqtZwpUyIoUeJXduyooWO6ebEZM9rRs+cyjh9vSvXqsRw7lmQ7krJAC5DKE88+u5xPP21FqVIx/PZbHcqUKWI7krLss8/a8uyzqzh5shG1asXq6TgvpAVI5bqXXlrJ1KkRlCy5nh076unoBirT5MmtefrpVZw40YSaNTdx5sx525FUHtICpHLVoEGrmTChBcWKbWL79jpafNSfTJvWhieeWMHRo82pVWuDzrDqRbQAqVwzcmQ077/flMDAWLZurUZwcGHbkZSH+vLLtjz00BIOH25J7dprSU1Ntx1J5QEtQCpXTJ68kWHD6hIQsIstWyprhwN1Tf/+d3v++tfF7N/fiiZNVuiICV5AC5DKcbNn76Bv38oUKnSEtWuDqVxZu1qr7Pnpp0huuWUxW7e2o23bJbbjqFymBUjlqF9+iaNLl+I4HEksWlSIevX0JlN1fVasaE/t2ktZuTKSzp0X246jcpEWIJVjtmw5SqdOIOLDd98l07p1RduRVD7kcBh+/bU1FSuu4qefIunbd4XtSCqXaAFSOSIhIZGWLRNISyvNtGnx3H13NduRVD7mdPqwZUtTihXbwKRJzRk7dr3tSCoXXLMAGWNeM8asNcacMcYcNcb8aIypn43tGhhjlhhjko0xB40xw40xOsNYAZSSkkaDBltITKzNkCGb6NXrmm8Ppa4pKMiPX3+tgp9fHAMGVGX2bJ3eu6DJTgsoEpgEtAJuA9KA/xljSl5pA2NMELAAOAI0B14EBgL9bzKv8jAul9C06UoSElrwyCMrGDWqpe1IqgAJDS3GkiVFcDiS6NKlKGvWHLIdSeWgaxYgEekkIp+LyGYR2QR0A0oDra+y2WNAYaB7xnbfAO8C/bUVVLDcc88SYmPbccsti/n663a246gCqGXL8syadRaXqzCRkWc5fPic7Ugqh9zINaCiGdudvMo6EcAyEck6zO3PQHkg7Ab2qTzQa6+tYc6cdlSsuIply7T4qNzz4IM1GTlyJykp1WnUaKveqFpA3EgBGg9sAFZdZZ0Q3KffsjqSZdlFjDHPGGOijTHRR48evYFIKq/NmBHL3/9enyJFYlm/vhG+vtqfReWuoUPDeeih5SQktCAiYpntOCoHXNenhjFmHNAGeEBErvUV5NLbmM0VnkdEpolIuIiEly6t9414upiYwzzxRHF8fE6xcmVpHWJH5Zl//7s9DRsuISYmksceW2o7jrpJ2S5AxpgPgK7AbSKy+xqrH+bPLZ0yGX9e2jJS+cipUym0a3ccl6soM2cm0rBhmWtvpFQOWrOmNaVLr+WrryKYOPFX23HUTchWATLGjAcexV18tmVjk1VAW2OMf5bnOgKHgL3XG1J5BpdLCA+PJjGxHoMHb+HBB2vajqS8kL+/L9HRNXE69/PSS+VZseKA7UjqBmXnPqCPgJ64Wz8njTEhGY/ALOuMNsYszLLZV0AS8IUxpr4x5n5gMDBORHSEwXzqgQeWsmtXGyIjFzN6tHa3VvZUrlyM774DEScdO57VGVXzqey0gPrg7vm2EIjP8hiQZZ1yQOat7yJyGneLpzwQDXwEjAXG5Uhqlefeey+G775rTblyq1mwQHu8KfvuvLMqb7yxg+TkWjRrtkFHz86HjKc1SMLDwyU6Otp2DJXFihUHaNu2ME7ncXbvDtGpFZRHufPOxcybF8lf/7qYn36KtB1HXcIYs05Ewi+3TPvOqqs6dSqFv/zlDCI+zJnj1OKjPM6cOe2pUGEVc+a0Yfz4DbbjqOugBUhdVUREFElJdRkyZBu33x5qO45Sf+JwGKKi6uN07qN///LExBy2HUllkxYgdUU9ey5j27Z2REQs1jHelEcrX74o//d/rozhehI4dy7VdiSVDVqA1GXNmrWdL75oTsmS61i8uK3tOEpd0z33VOfFFzdy9mxD2ra92kAtylNoAVJ/cuDAGbp188PhOMny5aE4nT62IymVLePHt6Jx4yVs2NCeV19dbTuOugYtQOoiLpfQqtVmLlyoyMSJR6lTJ9h2JKWuy7JlERQpsoX33qvD0qX7bcdRV6EFSF2kW7dl7N/fik6dltOnT0PbcZS6boGBTubODQLgrrvO6PUgD6YFSGWaNWs7X33VguDgtfz0k95sqvKvdu0qMWhQLImJ9WjXTq8HeSotQAqAhIREunVz4nCcZOnSKjq9gsr33n33Fho1WsL69e0ZPjzKdhx1GfopowBo2zaGCxdCGTv2iF73UQXG0qW3EBCwjZEjq7Jhgw7E72m0ACn691/Fjh1tadVqKS+/3Nh2HKVyTFCQH7NnOxEpzO23HyAtzWU7kspCC5CXW7HiAB98UJciRTazYEFr23GUynF33VWVbt2iOXGiGffeq5PYeRItQF4sJSWNu+46Dhh++imIwoUL2Y6kVK744ou2VKiwirlzW/PPf261HUdl0ALkxe65ZzlnzjSid++NREZWth1HqVzjcBiWLKmNj89RevUK0PmDPIQWIC81Y0YsCxa0plKllUyapKfeVMFXrVoJRo8+zIULVbjttrW24yi0AHmlEyeSefJJJw7HMRYtqoPDYWxHUipPDBzYlKZNF7NpU3tGjdJ5x2zTAuSFOnSIIjW1GqNGHaJatRK24yiVpxYuvAU/v98YPrwiO3eesB3Hq2kB8jLvvx/D+vXtadhwCYMHN7MdR6k8V7y4P59/nobLVZLbbtuuU3lbpAXIixw4cIbXXiuL07mbhQub246jlDVdu9bmL39ZyYEDEfTrp0P12KIFyIt06LCB9PQQPvooieDgwrbjKGXV99+3oUiRzUycWIvNm4/ajuOVtAB5iVGjotm+vR0tWy6jV6/6tuMoZZ2/vy8zZwYgUoS//GWXnoqzQAuQF4iLO83w4RXw8/uNefNusR1HKY9x993V+OtfVxMffwsvvLDSdhyvowXIC3TsuBGXqzRTp6ZSvLi/7ThKeZRvv21LYOAmJk+uqwOW5jEtQAXcqFHR7NzZllatltO9e13bcZTyOE6nD7NnByLiz1137dFTcXlIC1ABdujQWUaMCMHp3MV//xthO45SHqtTpyrcddca4uNvYeDA1bbjeA0tQAXYHXfEkJ5engkTkggK8rMdRymP9s03bShceCsffFBdb1DNI1qACqhJkzayaVNbGjVaRu/eDWzHUcrj+fv78vnnPogU5447dMTsvKAFqAA6dSqFfv0C8fE5xLx5OtqBUtn18MO1aNNmBbt3t2HkSB0rLrdpASqAOndeTWpqVd5++wghIYG24yiVr8yZE4HTuYs33gghISHRdpwCTQtQAfPDD7+xfHkrqlZdzmuvaetHqesVFOTHBx8kkp5ekbvu0mkbcpMWoAIkLc3F448nYcw5fvqptu04SuVbffo0pHbtpaxb14avv95mO06Bla0CZIxpZ4z5wRhz0Bgjxpge11g/LGO9Sx935EhqdVlPPrmcs2cb0qPHFurUCbYdR6l8be7chjgcJ+jVy0VqarrtOAVSdltAgcBm4CUg+Tpe/w6gXJbHL9eVTmXbxo0JTJ/ekOLF1/PJJ21sx1Eq36tSpTh9+/5GUlJdunZdbjtOgZStAiQic0VkiIjMBlzX8frHReRwlkfqjcVU13LPPTsRCeCrr4rpDKdK5ZAPP4wgOHgt337blLVr423HKXBy+xrQt8aYBGPMCmPMg7m8L6/13nsxxMW1pn37Vdx5Z1XbcZQqMBwOwzffhAC+3HffXttxCpzcKkDngAHAw8BdwEJgljHm8cutbIx5xhgTbYyJPnpU5+W4HmfOnOf110vg6xvHd9/pSNdK5bR27SrRocMaDh6MYNQovTcoJ+VKARKRYyIyVkRWi0i0iAwHpgKDrrD+NBEJF5Hw0qVL50akAuvBB1dx4UIVRow4qiNdK5VLvvkmgkKF9vDmm8GcOpViO06BkZfdsNcANfJwfwXe0qX7WbCgJRUqrOL118Ntx1GqwAoK8mPkyJNcuBDG/ffrYKU5JS8LUGNAr+LloIcfPgSk8803obajKFXgDRrUlEqVVrJo0S388kuc7TgFQnbvAwo0xjQ2xjTO2KZyxs+VM5aPNsYszLJ+d2PMo8aYOsaYWsaYAUBfYGJu/CO80fDhURw50pK77oqmZcvytuMo5RW++64qkErXrjpxXU7IbgsoHFif8QgA3sz4+1sZy8sB1S7Z5nUgGlgLPAI8KSIf3Gxg5R5sdPTosjidu/i//2ttO45SXqNp0xDuvjuGhIQWvPGGDtNzs4yIZ83+Fx4eLtHR2tPkav7yl8UsWBDJu+/GMGhQU9txlPIq586lUqrUfkR8OHasnM61dQ3GmHUictmL1DoWXD6zZs0hFixoTvnyq7X4KGVBYKCTN988xYULYTz00CrbcfI1LUD5zEMP7QUczJpV0XYUpbzW4MHNKFduNfPnN9cREm6CFqB8ZPz4Dezf34r27dfQpo0WIKVsmjmzAuDgwQf32I6Sb2kByidSU9MZPDgAH58DfPttS9txlPJ67dpVol27Nezb14p//ONX23HyJS1A+cSTT64gJaUWL764n5IlA2zHUUoB33zTAh+fQwwa5CQt7XrGaVagBShf2LfvNF99VYegoF8ZM0bHe1PKUwQHF6Z37z0kJ9fhmWdW2I6T72gBygfuv389IqWYNMmpUy0o5WEmTmxFYOAmvvyyJocOnbUdJ1/RAuThFizYy7p1rahRYwWPPVbHdhyl1CUcDsOECQ5crrI8+OA623HyFS1AHu6JJ44AKcyeXct2FKXUFfTsWY+qVZezatUtLF2633acfEMLkAcbO3Y9hw+3pFOnGBo2LGM7jlLqKmbNqgak8/jjWoCySwuQh0pLczFsmD8+PgeYOVM7Hijl6cLDy9G+/Vr272/FtGmbbMfJF7QAeag+fVaSnFyH3r3jdKI5pfKJmTPDcTiO8Mor4HJ51jibnkgLkAc6diyJTz+tRpEiWxg/PsJ2HKVUNoWEBNKt2w7OnWvAgAE6cd21aAHyQF27RuFylWP06DR8ffW/SKn8ZNq0Vvj57WTixPKcOXPedhyPpp9uHmbjxgT+9z/3QIcvvNDIdhyl1HVyOn0YNuw0aWmhdOumo2VfjRYgD9O16zbAn88/L2s7ilLqBg0dGk6pUtH8+GNj9uw5ZTuOx9IC5EF+/nkPW7e2okGDlXTqVMV2HKXUTZg0qSgiQTzyyAbbUTyWFiAP8tRTh4Fkvvqqru0oSqmb9PDDtahadSVRUbewZs0h23E8khYgD/Hpp5s5eDCCyMh11K9f2nYcpVQO+Oc/wwBDt267bEfxSFqAPIDLJbzyShrGHGXGjGa24yilckjr1hVp1mwVO3e24j//2Wk7jsfRAuQB3normtOnG/Pww1spX76o7ThKqRz09dcNgXM8++xJ21E8jhYgy1JT0/n734tTqNBePvusle04SqkcVqNGSTp1Wk9CQgudOfUSWoAse+ml1Zw/X4PnnjtE4cKFbMdRSuWC6dNb4HDEM2SI0SF6stACZNG5c6l88klFAgK2MXasDjiqVEEVHFyYLl12cvZsQ95+O9p2HI+hBciip59eTVpaKK++elaH3FGqgJs27RZ8feP4+9+DSEtz2Y7jEfRTz5Jjx5L4979rEhT0K8OGhduOo5TKZYGBTp5++gApKbV45RUdqBS0AFnzxBNRuFwhjBrlntJXKVXwTZgQgb//DiZPLkdS0gXbcazTAmTBvn2nmTevEaVLr+X553XAUaW8ha+vg/79T3LhQhWefVYHKtUCZMHjj69HpATjx+s9P0p5m7ffbkFg4Ca++qo6p06l2I5jlRagPLZz5wmWLWtKhQqr6Nq1tu04Sqk85nAYhg27QHp6eXr1WmM7jlVagPJYt24bgUA++kjHe1PKWw0Y0IRixTbwn//U5tixJNtxrMlWATLGtDPG/GCMOWiMEWNMj2xs08AYs8QYk5yx3XBjjFdfbd+8+Shr1jQnNHQV995b3XYcpZQlDodh5EiDy1WWnj2jbMexJrstoEBgM/ASkHytlY0xQcAC4AjQHHgRGAj0v7GYBUP37lsAf6ZOLW87ilLKsuefb0TJkuuYM6c+hw+fsx3HimwVIBGZKyJDRGQ2kJ07qB4DCgPdRWSziHwDvAv099ZWUHR0PDExt1C9uk42p5Rye+89P0SCeeIJ7xwdIbeuAUUAy0Qka2vpZ6A8EJZL+/RoPXvuAHz49NMw21GUUh7iqafqU6ZMFP/7XyP27TttO06ey60CFIL79FtWR7Is8ypr1hxi8+ZbqF17Fe3aVbIdRynlQcaNK4pICXr2XG87Sp7LzV5wlw75aq7wPMaYZ4wx0caY6KNHj+ZiJDt69doJOPjkk6q2oyilPMxjj9UhJGQNixY19rpWUG4VoMP8uaVTJuPPS1tGiMg0EQkXkfDSpQtW9+S1a+MzWz+tW1e0HUcp5YHGjAlCpDhPPuldraDcKkCrgLbGGP8sz3UEDgF7c2mfHumpp3agrR+l1NX83gr65RfvuhaU3fuAAo0xjY0xjTO2qZzxc+WM5aONMQuzbPIVkAR8YYypb4y5HxgMjBMRr5mNae3aeDZt0taPUura3n8/CJESXtUKym4LKBxYn/EIAN7M+PtbGcvLAdV+X1lETuNu8ZQHooGPgLHAuBxJnU/83vr5+GNt/Silru7xx+tQtqy7FXTgwBnbcfJEdu8DWiwi5jKPHhnLe4hI2CXbbBKRdiLiLyLlRORNb2r9REf/0fpp00ZbP0qpa3NfCypBz54xtqPkCR0LLpc8/bT7vp+pU/WmU6VU9jz+eB3KlIli4cIGXjE6ghagXLB581E2bGhB9ep6349S6vq8805hRErx1FMFf3QELUC54KmntgBOJk3S4qOUuj5PPVWfkiXX8d//1uPEiWsOvZmvaQHKYTt3niAqKpzQ0NV07BhmO45SKh96801fRErz9NMFe6RsLUA57Mkn3fP9jB9f1nYUpVQ+9fzzjShWbAPffVeTM2fO246Ta7QA5aC4uNMsX96EChV0vh+l1M0ZMsSFy1WO554ruLOmagHKQU89tR4oxtixJW1HUUrlcwMGNCEwcBOzZlUlKemC7Ti5QgtQDklISOSXXxpQpkwUXbrUsh1HKZXPORyGfv2SSU+vyMsvF8xWkBagHPLss9GIlOKttwJsR1FKFRBvvNEcf//tfPllOdLSsjMXaP6iBSgHnDuXyg8/1KBYsQ307t3AdhylVAHhcBh69TpOamo1hg4teD3itADlgBdeWEN6enkGDkyzHUUpVcC8/34LfH3j+OijorhcBWs0My1ANyk1NZ0ZMyoSEBDLa681sx1HKVXA+Pv70rVrHImJ9RgzpmCNlK0F6CYNHhzFhQtVeO65Uzgc5tobKKXUdfrHP1ricMTz97/bTpKztADdBJdLmDKlOE7nbkaPbmE7jlKqgAoK8uPuuw3LSjoAABv6SURBVLdz8mRTPvlks+04OUYL0E0YPXodycl1ePTRgzidPrbjKKUKsKlTwzHmJMOHF5xRsrUA3YSxYx04HPGMH6+tH6VU7goJCaR161+Jj2/B3Lm7bcfJEVqAbtA//7mVkyebcued2wkK8rMdRynlBaZMqQecp1+/g7aj5AgtQDdoyJDTwGmmTGlqO4pSykvUq1ea+vWj2LGjBTExh23HuWlagG7AwoVxHDzYkoiI9VSsGGQ7jlLKi0ycWAXw5bnnttmOctO0AN2AF16IA1KZMqWu7ShKKS8TGVmZSpXWEBXVhH37TtuOc1O0AF2njRsTiI1tQe3aUTRsWMZ2HKWUF3r33ZJAMXr3zt83pmoBuk59+mwFnEyYUNl2FKWUl+ratTYlS65j/vw6+XrCOi1A1yEhIZGVKxtRvnyUTretlLJq4EBwucrSv/9a21FumBag69C3bzQiJRgxoojtKEopLzdoUFP8/bczfXrZfDtIqRagbEpNTee776oQGLiJZ57RKReUUnY5HIZu3Y5y/nwN3nlnne04N0QLUDa99loUaWmVee65RNtRlFIKgHHjWuBwxDNuXP78KM+fqfOYyyVMnVoUX984Ro5sbjuOUkoBEBjopFMn9yCls2Zttx3numkByoapUzeRmFifBx7Yq4OOKqU8yuTJTYCzvPbaUdtRrpsWoGx4++0kjDnOP/6hrR+llGcJDS1G06br2LOnJWvWHLId57poAbqGhQvjiI9vQevWmwgOLmw7jlJK/cmECdUBw4sv7rAd5bpoAbqGfv32AmlMnFjHdhSllLqs1q0rUrHiWqKimnD4cP6ZLyjbBcgY08cYs8cYk2KMWWeMaXuVdcOMMXKZxx05EztvxMWdZtOmZlSrFkXjxmVtx1FKqSsaMaIoUIwXXsg/XbKzVYCMMV2A8cA7QBNgJfBfY8y1xqO5AyiX5fHLjUfNe336rAcCeeed0rajKKXUVfXqVZ/AwE18/30oqanptuNkS3ZbQP2BL0TkYxGJFZEXgHjguWtsd1xEDmd5pN5U2jyUlHSBn3+uSfHi63n44Vq24yil1DX16nWWCxfCeOONaNtRsuWaBcgY4wSaAfMvWTQfaHWNzb81xiQYY1YYYx68wYxWDB68lvT08rzwwgXbUZRSKltGj26Bj88BJk3KH7M0Z6cFFAz4AEcuef4IEHKFbc4BA4CHgbuAhcAsY8zjN5gzT7lcwmefFaNQoT0MHx5uO45SSmWLv78vd975G6dPN2bGjFjbca7J9zrWvXS0O3OZ59wrihwDxmZ5KtoYEwwMAqZfur4x5hngGYDKle1Pc/Dxx5tJTGxAly5L8PWtYjtOvudyuTh27BinTp0iPT1/nJtWecvHx4fixYsTHByMw6Gdc2/GRx814aefzjJs2HEee8x2mqvLTgE6BqTz59ZOGf7cKrqaNUDPyy0QkWnANIDw8HDrw7qOGnUWY04xYYK2fnLCgQMHMMYQFhZGoUKFMMbYjqQ8iIhw4cIFjhw5woEDBzziS2h+VrlyMRo2XMLGjRFs3Jjg0RNnXvOrRkbHgXVAx0sWdcTdGy67GuPuuODR1qw5xP79LWjWbANlyui0CzkhMTGRChUq4HQ6tfioPzHG4HQ6qVChAomJOthvThgzJhTw5cUXt9qOclXZPQU3DviXMSYKWAE8C5QHpgAYY0YDLUTk9oyfuwMXgPWAC+gM9AVezdH0ueDll3cAZfngg+q2oxQoelpFXYu+R3JOx45hlCkTxdKl9Thz5jxBQZ7ZKSFb/+MiMgt4GXgd2AC0Ae4SkbiMVcoB1S7Z7HUgGlgLPAI8KSIf5ETo3HLiRDJr1jSkXLm1tGlT0XYcpZS6YS+/7INIaQYN8twZU7P9lUNEJolImIj4iUgzEVmaZVkPEQnL8vOXIlJXRIqISJCIhIvInzofeJp+/aIRKcmrr/rbjqKUUjfl1Veb4nTu4l//KumxM6ZqmzeDyyXMmhWCv/92Xnihke04SgEQFhbGmDFjbMdQ+ZDDYXjggYMkJdVl6tRNtuNclhagDB98sIHz52vQtWsCDodeKFdK5X8fftgMY07xzjue2blDC1CGMWNSMeYY48bpnD/qDyLCe++9R7Vq1QgICKBBgwZMn+4+m/zUU09Rr149kpOTAUhPT6dNmzbcfffdmdsPHjyYWrVqERAQQFhYGIMGDSIlJeWifcyZM4eWLVsSEBBAqVKl6Ny5MykpKURGRhIXF8fAgQMxxmgPQnXdypQpQrNmGzhwoLlHzhWkBQhYvvwAhw+HExGxmeLF9fqP+sPrr7/Op59+ykcffcTWrVt57bXX6N27N3PmzGHChAlcuHCBAQMGADBq1Ch+++03Pvvss8ztixQpwmeffUZsbCyTJk1i5syZjBo1KnP5vHnzuPfee+nYsSPr1q1j0aJFtG/fHpfLxbfffkvFihUZPnw48fHxxMd7/F0MygO5e/Q6eOUVz5sryIh41sWp8PBwiY7O24H0WrZcTFRUG1avTqBly/J5um9vEBsbS506F8+n9PLLL7Nhw4Y8zdG4cWM+/PDDbK+fmJhIcHAw8+fPp23bP2Yfefnll9mxYwdz584lOjqaVq1aMXjwYEaPHs0PP/zAnXfeecXXnDJlCmPGjOG3334DoHXr1lSqVImZM2dedv2wsDCef/75zCJX0F3uvaJuXrlyazhypConThTN8y/Zxph1InLZu/q9vgV04kQya9c2pEKFtVp81EW2bt1KSkoKd9xxB4GBgZmPyZMns2vXLgDCw8MZOnQob7/9Ns8888yfis/s2bNp06YNISEhBAYG0q9fP/bt25e5fP369dx+++15+u9S3qdfv0KIlGbwYM8aJft6xoIrkAYMiEakLQMGBNiO4lWupyVii8vlAuDHH3/80/AwhQoVAtzXiJYvX46Pjw+7du1CRDKv1axevZpHHnmEESNG8MEHH1C8eHF++OEHr2nNKM8xYEAThg3bzfTpxZkyxXaaP3h1C8jlEmbODMbPbycvvqhdr9XF6tati5+fH3FxcVSvXv2iR2hoKADjxo0jJiaGpUuXsnr1aiZOnJi5/YoVK6hQoQLDhg2jefPm1KhRg7i4uIv20aRJExYuXHjFDE6nUwdwVTfN4TDcc89+EhPr8/nnW2zHyeTVLaCPP95McnIDHn10KQ5HDdtxlIcpWrQoAwYMYMCAAYgI7dq149y5c6xevRqHw0HLli0ZOnQoM2bMoFWrVkyePJknn3yS2267jfr161OzZk0OHjzIjBkziIiI4Oeff+brr7++aB9Dhw6lc+fOVK9enUcffRQRYf78+fTu3ZvChQsTFhbGsmXLePzxx/Hz8yM4ONjS0VD53fjxTZk9+wwjR56k52WHhbZARDzq0axZM8kroaHLBU5JfPzZPNunN9q6davtCDfM5XLJhAkTpE6dOuJ0OiU4OFg6dOggP/zwg9SrV0969ux50fqPP/64NGzYUFJSUkREZPDgwRIcHCxFihSR++67TyZNmiTuX7s/fP/999K0aVNxOp1SqlQp6dy5syQnJ4uIyKpVq6Rhw4bi5+f3p+0Kovz8XskPGjRYLJAimzcn5Nk+gWi5wue91/aC27gxgUaNStCkyUpiYtrn+v68mfZsUtml75XcNXfubv7616p06LCYBQsi82Sf2gvuMvr33woU4t13Q21HUUqpPHHXXVUpWXIdixbVIDXV/rVFryxAKSlpLF5ck1KlounYMcx2HKWUyjO9el0gPb0Cb75pv0u2VxagN99cR3p6eZ55xv43AKWUyktvvhmOwxHP1Kk+tqN4ZwGaNs0XH5+DDB/ezHYUpZTKU/7+vkRGbuf48XAWLoy79ga5yOsK0M8/7+HEiWbceutO/P29uhe6UspLjR1bB7jAq6/usZrD6wrQkCFxQCpjx9a1HUUppaxo3LgslSqtJSamESdOJFvL4VUF6NixJNavb0zlytE0bFjGdhyllLLmlVcKI1KCQYPWWcvgVQVo4MB1iBRn4MAitqMopZRVL7zQCD+/3/j665LWMnhVAfr3v0vh57eTPn0a2o6ilFJWORyGzp3dU3ZPnx5rJ4OVvVowfXosSUl16dz5kE65rbzCF198QWBgoO0YyoONHdsYSGTUqKNW9u81Bch9gBMzDrhSBV+XLl3YvXt3ru9n8eLFGGM4duxYru9L5azKlYtRs2YM27Y1Zd++03m+f68oQPv2nWbbtqbUrBlD5crFbMdR6qpcLleOTMEQEBBAmTLa2UZd3dChpYBABg3K2xmKwUsKkPvABmYcaKWyb968ebRt25YSJUpQsmRJOnXqRGys+3z53r17Mcbw1Vdf0aZNG/z9/alduzbz58/P3P731sFPP/1E48aN8ff3p1mzZqxb90fPo99Plc2dO5f69evjdDqJjY3l5MmTdO/enRIlShAQEECHDh3YssU9l8vRo0cpV64cb731VubrbNy4EX9/f2bPnn3R6/7ujTfeoH79+nz55ZeEhYURGBhIz549SU1NZdKkSVSqVIlSpUrRv3//zMn4AKZPn07z5s0pWrQoZcqU4aGHHuLgwYOZx+DWW28FoHTp0hhj6NGjB+Aeaf+9996jWrVqBAQE0KBBA6ZPn56T/z0qBzzxRF0KF97K99+Xw+XK48GprzRMtq1HTk/HkJ7uEn//7RIQoMO825Kfh9ifPXu2zJ49W3bs2CG//vqrPPTQQ1KtWjU5f/687NmzRwCpUKGCzJo1S2JjY+X5558Xf39/OXDggIiILFq0SACpVauWzJs3TzZt2iQPPviglC1bVhITE0VE5PPPPxcfHx+JiIiQ5cuXy/bt2+XMmTNyzz33SK1atWTJkiWyceNG6dy5s1SsWFGSkpJERGTevHlSqFAhWblypSQlJUndunWlR48emdk///xzKVKkSObPI0aMyJwWYtOmTTJv3jwpUqSI3HHHHdKjRw/ZunWrfPvtt+Lr6yuzZ8/O3O7TTz+VOXPmyK5du2TNmjUSGRkpbdu2FRGRtLQ0+eabbwSQLVu2SHx8vJw6dUpERIYMGSI1a9aU//73v7J7926ZMWOGFC5cWH766acrHu/8/F7Jz7p1WyIgMmXKxhx/ba4yHYP1gnPpI6cL0JQpGwXcB1jZcbkPlZdeEmnfPm8fL7108/+Wc+fOicPhkGXLlmUWoJEjR2YuT09Plxo1asjQoUNF5I8CNH369Mx1zp49K8WKFZOPP/5YRNyFApDo6OjMdXbs2CGALFnyx/v21KlTEhQUlLmdiMhLL70kVapUkR49eki1atXk7Nk/5ra6XAHy9/fPLBAiIg888IAEBwfL+fPnM59r37699O3b94rHIDY2VgDZv3//Rf/Go0ePXnSc/P39ZenSpRdt+9JLL8mdd955xdfWAmRHfPxZgdNSteqyHH/tqxWgAn8K7r33TgNneO+9prajqHxo165dPProo1SrVo2goCDKli2Ly+Vi3759metERERk/v33mVK3bt160etkXScwMJAGDRpctI6vry+NG//RQSY2NhaHw3HRdsWKFfvTdu+++y5Op5N//vOfzJgx45q93ipXrkyxYn9cBy1btiw1a9bE6XRe9FxCQkLmzzExMdx7772EhoZStGhRwsPdU7tkPQaX2rp1KykpKdxxxx0EBgZmPiZPnsyuXbuumlHlvZCQQOrXX8/u3eHs3Hkiz/ZboAdD27XrJLt3h1OvXhQhIe1sx1FZfPih7QTZ07lzZypUqMDUqVOpUKECvr6+1K1bl9TU1Bzdj5+fHz4+f4xOLFeZKNKYP24j2Lt3L/v378cYw+7du2nZsuVV91OoUKE/vdblnvu9E0RiYiKdOnWiQ4cO/Otf/6JMmTIcO3aMtm3bXvUY/H4N6ccff6Ry5cpXzaA8w4gR5XjoIX8GDFjN999H5sk+C3QLaODAjYA/w4eXtR1F5UPHjx8nNjaWIUOG0KFDB+rUqcPZs2dJS0u7aL3Vq1dn/l1EiIqK+tOsnlnXSUxMZPPmzVed+bNu3bq4XC5WrVqV+dyZM2fYtGkTdeu6xzG8cOECjz32GPfccw9jxozhueeeu2qr5EZs27aNY8eO8c4779CuXTtq1659UesIyGw9Ze25V7duXfz8/IiLi6N69eoXPUJDdRJIT/TggzUJDNzEvHmV8qwzQoFtAblcwty5FShSZDMPP1zfdhyVD5UoUYLg4GA+/vhjKlWqxMGDBxk4cCC+vhf/2kyePJmaNWvSoEEDJk2aRFxcHM8999xF64wcOZLSpUtTvnx53nrrLZxOJ48++ugV912jRg3uvfdeevfuzbRp0yhevDhDhw4lKCgoc7thw4aRkJDA//73P4oVK8a8efPo1q0bixYtwuHIme+WlStXxs/Pj3/84x/07duX2NhYhg0bdtE6oaGhGGOYM2cOnTt3JiAggKJFizJgwAAGDBiAiNCuXTvOnTvH6tWrcTgcPPPMMzmST+Wshx46xeefN2DSpF95/vlGub/DK10csvXIqU4Iv3c+6N596bVXVrkqP19YXrhwodSrV0/8/PykXr16mT3HPv/888xOCNOnT5eIiAjx8/OTmjVryty5czO3//0C/ffffy8NGjQQp9MpTZo0kaioqMx1Lu0s8LsTJ07IE088IcWLFxd/f3+5/fbbZfPmzSIisnjxYvH19ZVFixZlrh8fHy+lS5eWd95557KvO2LECKlXr95F++jbt6+0b9/+oue6dOkiDzzwQObPM2fOlKpVq4qfn580b95c5s2bJ8BF+37rrbckJCREjDHSvXt3ERFxuVwyYcIEqVOnjjidTgkODpYOHTrI/Pnzr3i88/N7pSA4cuScwKkc7YzAVTohGLnKuWYbwsPDJTr65qeKrV59Obt2NSI+3hASosOR2BQbG3vV00351d69e6lSpQpr167NvDB/qcWLF3Prrbdy9OhRgoOD8zhh/lNQ3yv5Sf36S9mypQW//ZZMtWolbvr1jDHrROSyvyDZbqcbY/oYY/YYY1KMMeuMMW2vsX4DY8wSY0yyMeagMWa4yXr1NBfFxZ1m166m1KmzXouPUkpdh2HDygL+DBq0Mdf3la0CZIzpAowH3gGaACuB/xpjKl9h/SBgAXAEaA68CAwE+udA5msaOHADUJihQ0vnxe6UUqrA6NKlFkWKbGHOnPK53hkhuy2g/sAXIvKxiMSKyAtAPPDcFdZ/DCgMdBeRzSLyDfAu0D+3W0Eul/DjjyEULryVxx7TprzKPWFhYYjIFU+/AURGRiIievpN5SsPPHCC8+dr8Mknm3N1P9csQMYYJ9AMmH/JovlAqytsFgEsE5Gsc73+DJQHwq4/ZvZ9+eVWUlJqcd99OjKvUkrdiHffbQKcZcyYU7m6n+y0gIIBH9yn07I6AoRcYZuQK6z/+7KLGGOeMcZEG2Oijx69uXkpdu9OxOncxXvvNbmp11FKKW8VEhJIw4Yxub6f67lZ4NKTgeYyz11r/cs9j4hME5FwEQkvXfrmrtu8/XYLkpOrUr580Zt6HZWzPK23pfI8+h7xLOvXt2PHjqv2Nbtp2SlAx4B0/txyKcOfWzm/O3yF9bnKNjlGZzz1LIUKFSI5OfnaKyqvlpycrMP0eJC8+By9ZgESkVRgHdDxkkUdcfeGu5xVQFtjjP8l6x8C9l5/TJWflSlThoMHD5KUlKTfctWfiAhJSUkcPHhQJ9DzMtkdimcc8C9jTBSwAngWd4eCKQDGmNFACxG5PWP9r4ARwBfGmJFATWAw8KboJ5DXCQoKAuDQoUNcuHDBchrliQoVKkTZsmUz3yvKO2SrAInILGNMKeB1oBywGbhLROIyVikHVMuy/mljTEfgIyAaOAmMxV3IlBcKCgrSDxel1EWyPRipiEwCJl1hWY/LPLcJ0DkQlFJKXVaBno5BKaWU59ICpJRSygotQEoppazQAqSUUsoKj5sPyBhzFIi75opXF4z7Blqlx+J3ehz+oMfCTY/DH3LzWISKyGWHuPG4ApQTjDHRV5oAydvosXDT4/AHPRZuehz+YOtY6Ck4pZRSVmgBUkopZUVBLUDTbAfwIHos3PQ4/EGPhZsehz9YORYF8hqQUkopz1dQW0BKKaU8nBYgpZRSVhS4AmSM6WOM2WOMSTHGrDPG5O6Ufh7GGPOaMWatMeaMMeaoMeZHY0x927lsM8YMMcaIMeYftrPYYIwpZ4z5MuM9kWKM2WqMaW87V14zxvgYY97O8hmxxxgz0hiT7YGZ8yNjTDtjzA/GmIMZvwc9LllujDFvGGMOGWOSjTGLjTH1cjtXgSpAxpguwHjgHaAJ7gnz/muMqWw1WN6KxD1qeSvgNiAN+J8xpqTNUDYZY24BngY22s5igzGmOO55vAzwV6AO8AKQYDOXJa8CfYEXgdrASxk/v2YzVB4IxD2NzkvA5aYnHgS8gvt90Rz3e2OBMaZoboYqUJ0QjDFrgI0i8nSW53YCs0WkoL/BLssYEwicBv4mIj/azpPXjDHFgBjcBWg4sFlEnrebKm8ZY94B2otIa9tZbDPG/AQcF5HuWZ77EiglInfbS5Z3jDHngOdF5IuMnw3u2ar/ISKjMp4LwF2EBojI1NzKUmBaQMYYJ9AMmH/Jovm4WwPeqiju/+eTtoNYMg33F5BfbAex6G/AGmPMLGNMgjFmgzHm+YwPHm+zHLjVGFMbwBhTF/eZgrlWU9lVBQghy2eniCQDS8nlz86CdN4zGPABjlzy/BGgQ97H8RjjgQ3AKttB8pox5mmgOtDNdhbLqgJ9gA+AvwONgYkZy7ztmti7uL+UbTXGpOP+DByVMeGmtwrJ+PNyn50VcnPHBakA/e7Sc4rmMs95BWPMOKAN0EZE0m3nyUvGmFq4rwW2FZFU23kscwDRWU5DrzfG1MB97cPbClAX4AngUWAL7mI83hizR0Q+tZrMvjz/7Cwwp+Bwj+Sazh/V/Hdl+HNlL/CMMR8AXYHbRGS37TwWROBuFW82xqQZY9KA9kCfjJ/97MbLU/HA1kueiwW8qXPO794HxojITBHZJCL/AsZR8DshXM3hjD/z/LOzwBSgjG+564COlyzqiLs3nNcwxozH/Q3vNhHZZjuPJd8BDXB/w/39EQ3MzPi7N7WKVgC1LnmuJjc/7Ul+VBj3F9Ws0ilAn4U3YA/uIpT52WmM8QfaksufnQXtFNw44F/GmCjcv3TPAuWBKVZT5SFjzEe4r3n8DThpjPn9W805ETlnL1neEpFTwKmszxljEoETIrLZTiprPgBWGmOGArNw36LwIjDEaio7fgQGG2P24D4F1wToD/zTaqpcltEbtnrGjw6gsjGmMe7fh33GmA+BocaYbcAO4HXgHPBVruYqSN2wwX0jKu4+7eVw93vvJyJL7abKO8aYK/2Hvikib+RlFk9jjFmMF3bDBjDG/BX3NbFawD7c134mSkH7ALiGjPta3gbuw32KKR53q/gtEUmxmS03GWMigUWXWfSliPTI6BE5AugNlADWAH1z+8tagStASiml8gdvPu+plFLKIi1ASimlrNACpJRSygotQEoppazQAqSUUsoKLUBKKaWs0AKklFLKCi1ASimlrNACpJRSyor/B5k1xaI8nNwoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# First, let's check that when the drag coefficient is zero, the solution is the\n",
    "# same as with the analytic one\n",
    "\n",
    "r_t0 = (0, 0, 0)\n",
    "v_t0 = initial_velocity(10, 45)\n",
    "t_f = flight_time(10, 45)\n",
    "\n",
    "# Let's use 100 time steps, to guarantee numerical convergence\n",
    "ts = np.linspace(0, t_f, 100)\n",
    "\n",
    "# compute the positions using the exact, analytic formula\n",
    "xs_a, ys_a, zs_a = parabolic_positions(ts, r_t0, v_t0)\n",
    "# compute the positions using the approximate, numeric method\n",
    "xs_n, ys_n, zs_n = parabolic_positions_numeric_drag(ts, r_t0, v_t0, c_d=0)\n",
    "\n",
    "# plot\n",
    "plt.rcParams.update({'font.size': 14})\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.plot(xs_a, ys_a, 'k-', label='exact')\n",
    "ax.plot(xs_n, ys_n, 'b-', label='approximate')\n",
    "ax.legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-01-31T23:38:25.489504Z",
     "start_time": "2021-01-31T23:38:25.306492Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a5324fec138a4acab510935acfef0fb3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='drag', max=1.0), Output(layout=Layout(height='350px'…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Now, let's explore the impact of the drag coefficient\n",
    "# Let's see that with an interactie plot\n",
    "\n",
    "# common parameters that will not change\n",
    "r_t0 = (0, 0, 0)\n",
    "v_t0 = initial_velocity(10, 45)\n",
    "t_f = flight_time(10, 45)\n",
    "\n",
    "ts = np.linspace(0, t_f, 100)\n",
    "\n",
    "# parameters that we want to change\n",
    "min_drag = 0\n",
    "max_drag = 1\n",
    "change_drag = 0.1\n",
    "\n",
    "\n",
    "def f(drag):\n",
    "    # compute the positions using the exact, analytic formula\n",
    "    xs_a, ys_a, zs_a = parabolic_positions(ts, r_t0, v_t0)\n",
    "    # compute the positions using the approximate, numeric method\n",
    "    xs_n, ys_n, zs_n = parabolic_positions_numeric_drag(ts, r_t0, v_t0, drag)\n",
    "\n",
    "    # plot\n",
    "    plt.rcParams.update({'font.size': 14})\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(1, 1, 1)\n",
    "    ax.plot(xs_a, ys_a, 'k-', label='no-drag')\n",
    "    ax.plot(xs_n[ys_n >= 0], ys_n[ys_n >= 0], 'b-', label='drag')\n",
    "    ax.set_xlabel(r\"X [m]\"), ax.set_ylabel(r\"Y [m]\")\n",
    "    ax.legend(loc=1)\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "# Output interactive plot\n",
    "interactive_plot = interactive(f, drag=(min_drag, max_drag, change_drag))\n",
    "output = interactive_plot.children[-1]\n",
    "output.layout.height = '350px'\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See how air drag reduces both the reach of the ball, and its maximum height. The numeric approach allows us to solve this problem with little complexity added relative to the case without drag.\n",
    "\n",
    "Ideally, you would code separate functions for the integration and the evaluation of the acceleration, so that you can use the same integration function with different accelerations.\n",
    "\n",
    "In practice, you would not code your own integration function, but use one from a library, since there are implementations that have been thoroughly tested for numeric accuracy and computing time efficiency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}