{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 9 \n",
    "<!-- ToC -->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-While statements: [>>](#While-statements)  \n",
    "--Basic structure: [>>](#Basic-structure)  \n",
    "--Another example: [>>](#Another-example)  \n",
    "--Use while or for loop?: [>>](#Use-while-or-for-loop?)  \n",
    "--A more complicated while loop: [>>](#A-more-complicated-while-loop)  \n",
    "--Projectile with air resistance re-visited: [>>](#Projectile-with-air-resistance-re-visited)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## While statements\n",
    "\n",
    "The use of the `while` statement is illustrated in the examples below.\n",
    "\n",
    "### Basic structure\n",
    "\n",
    "This example shows the syntax of a `while` loop. (As we will see below, it can be written in a more condensed way, and a `while` loop is probably not the best way of tackling this problem, but it does show how `while` loops work!) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " \n",
      "After 11 loops, thingToCalc = 518918400.0\n"
     ]
    }
   ],
   "source": [
    "#\n",
    "# Use debug to turn debug printing on and off\n",
    "#debug = True\n",
    "debug = False\n",
    "#\n",
    "# Initialise variables\n",
    "Test = True\n",
    "nLoopMax = 10\n",
    "nLoop = 0\n",
    "thingToCalc = 13.0\n",
    "#\n",
    "# Define the while loop\n",
    "while Test:\n",
    "    #\n",
    "    # Calculate the quantity(ies) you are interested in\n",
    "    thingToCalc = thingToCalc + nLoop*thingToCalc\n",
    "    if debug:\n",
    "        print(\"nLoop\",nLoop,\"thingToCalc =\",thingToCalc)\n",
    "    #\n",
    "    # Update the variables used to determine whether the while loop should continue\n",
    "    Test = nLoop < nLoopMax\n",
    "    nLoop = nLoop + 1\n",
    "#\n",
    "# Use results\n",
    "print(\" \")\n",
    "print(\"After\",nLoop,\"loops, thingToCalc =\",thingToCalc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no reason to introduce an explicit `Test` variable, as is shown in the next example.\n",
    "\n",
    "### Another example\n",
    "\n",
    "Print a string vertically instead of horizontally. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Initialise variables\n",
    "myString = \"Constantinople\"\n",
    "nLoopMax = len(myString)\n",
    "nLoop = 0\n",
    "#\n",
    "# Define the while loop\n",
    "while nLoop < len(myString):\n",
    "    #\n",
    "    # Calculate the quantity(ies) you are interested in\n",
    "    print(f\"{myString[nLoop]}\")\n",
    "    #\n",
    "    # Update the variables used to determine whether the while loop should continue\n",
    "    nLoop = nLoop + 1\n",
    "#\n",
    "# Use results\n",
    "print(\" \")\n",
    "print(f\"Number of characters in '{myString}' is {nLoop}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use while or for loop?\n",
    "\n",
    "The `while` loop is usually not the best way of tackling a problem where you can calculate the number of iterations you need. The `for` loop is designed for this purpose and is often a bit quicker and is usually easier to understand!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C\n",
      "o\n",
      "n\n",
      "s\n",
      "t\n",
      "a\n",
      "n\n",
      "t\n",
      "i\n",
      "n\n",
      "o\n",
      "p\n",
      "l\n",
      "e\n",
      " \n",
      "Number of characters in 'Constantinople' is 14\n"
     ]
    }
   ],
   "source": [
    "#\n",
    "# Initialise variables\n",
    "myString = \"Constantinople\"\n",
    "nLoopMax = len(myString)\n",
    "#\n",
    "# Define the for loop\n",
    "for nLoop in range(0, nLoopMax):\n",
    "    #\n",
    "    # Calculate the quantity(ies) you are interested in    \n",
    "    print(f\"{myString[nLoop]}\")\n",
    "#\n",
    "# Use results\n",
    "print(\" \")\n",
    "print(f\"Number of characters in '{myString}' is {nLoopMax}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A more complicated while loop\n",
    "\n",
    "Calculate the sum of cubes of the integers, stopping when the sum reaches a chosen maximum value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " \n",
      "Print all cubes and their sum up to sum value of 3000\n",
      "n \t n**3 \t Sum n**3\n",
      "0 \t 0 \t 0.0\n",
      "1 \t 1 \t 1.0\n",
      "2 \t 8 \t 9.0\n",
      "3 \t 27 \t 36.0\n",
      "4 \t 64 \t 100.0\n",
      "5 \t 125 \t 225.0\n",
      "6 \t 216 \t 441.0\n",
      "7 \t 343 \t 784.0\n",
      "8 \t 512 \t 1296.0\n",
      "9 \t 729 \t 2025.0\n",
      "10 \t 1000 \t 3025.0\n",
      " \n",
      "Calculation of sum of cubes completed.\n"
     ]
    }
   ],
   "source": [
    "#debug = True\n",
    "debug = False\n",
    "# \n",
    "# Initialise variables\n",
    "nCubeSum = 0.0\n",
    "n = 0\n",
    "nCubeSumMax = 3000\n",
    "#\n",
    "#  Print columns headers\n",
    "print(\" \")\n",
    "print(\"Print all cubes and their sum up to sum value of\",nCubeSumMax)\n",
    "print(\"n \\t n**3 \\t Sum n**3\")\n",
    "# \n",
    "# Start the loop, which will run as long as the test in the while statement is True \n",
    "#\n",
    "if debug:\n",
    "    #\n",
    "    # Look at test variables that will be used in first iteration of loop\n",
    "    print(f\" \")\n",
    "    print(f\"Initial test, loop {n}\")\n",
    "    print(f\"nCubeSum is {nCubeSum} so nCubeSum < nCubeSumMax is {nCubeSum < nCubeSumMax}\")\n",
    "    print(\"n \\t n**3 \\t Sum n**3\")\n",
    "#\n",
    "while nCubeSum < nCubeMax:\n",
    "    #\n",
    "    # The loop is indicated by the indentation. Do the required calculations in the loop\n",
    "    nCubeSum = nCubeSum + n**3\n",
    "    print(f\"{n} \\t {n**3} \\t {nCubeSum}\")\n",
    "    # \n",
    "    # Update the value of the variable used in the while test\n",
    "    n = n + 1\n",
    "    if debug:\n",
    "        #\n",
    "        # Look at test variables that will be used in next iteration of loop\n",
    "        print(f\" \")\n",
    "        print(f\"Test for loop {n}\")\n",
    "        print(f\"nCubeSum is {nCubeSum} so nCubeSum < nCubeSumMax is {nCubeSum < nCubeSumMax}\")\n",
    "        print(\"n \\t n**3 \\t Sum n**3\")\n",
    "    #\n",
    "#\n",
    "# The loop finishes where the indentation stops\n",
    "print(\" \")\n",
    "print(\"Calculation of sum of cubes completed.\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Projectile with air resistance re-visited\n",
    "\n",
    "Calculate the path of a projectile with air resistance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " \n",
      "Radius 0.025 m, mass 0.131 kg.\n",
      "Initial position ( 0.00,  0.00) m.\n",
      "Initial velocity (30.00, 40.00) m/s.\n",
      "Maximum time 8.155 s, time step 8.155e-03 s.\n",
      " \n",
      "Final step is number 836, actual flight time 6.818 s.\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAFNCAYAAAAJsbjVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyZklEQVR4nO3dedzVY/7H8denfWNSkZSpEJLJcjNiLClZosJvJltkzZYtxmQwjD1LlrFFUmPfBskuZMlWaKamUhKVEBVRtNyf3x/Xt3HLfXev51zne877+Xicx9m+53ve9133+ZzvdV3f6zJ3R0RECk+t2AFERCQOFQARkQKlAiAiUqBUAERECpQKgIhIgVIBEBEpUCoAkhfM7Hsz2yST+zWzkWZ2WU2/R3WZ2W5mNj12DkkfFQCJxsxmm9my5EP2SzO728yaVGVf7t7E3WdVM8+rZnZ8Te+3ppmZm9lmq++7++vuvkXMTJJOKgASWy93bwJsD+wIXLDmBmZWJ+upRAqACoDkBHefBzwLbA3/+5Z7qpnNAGYkj51gZjPNbKGZjTazjVa/vuS3YjOrb2bXmtlnyZHF7WbWsMS2fczsQzP7zsw+NrN9zexyYDfg5uSI5OY197smMzsg2c9iMxtvZp3L+vmS/ZxiZjPMbImZXWpmm5rZW0mOh82sXontS/1Zzey1ZJNJSc5DzKyrmc0t8dqOydHMYjObYma9Szw30sxuMbOnkxzvmNmmFf13kjzj7rroEuUCzAb2Sm5vDEwBLk3uO/Ai0AxoCHQDviYcKdQH/gG8VmJfDmyW3L4BGJ28dh3gKeDK5LnfA98CPQhfgFoDWybPvQocv0bGkvsdCVyW3N4e+ArYCagN9E9+nvpl/KyeZFoX6AT8BIwFNgF+A/wX6J9sW+GfNbnfFZib3K4LzAT+CtRL9rUE2KLEz7Aw+T3UAe4DHoz9f0GXOBcdAUhsT5jZYuANYBxwRYnnrnT3he6+DDgCGOHu77v7T8B5wM5m1q7kzszMgBOAs5LXLkn2eWiyyXHJfl5092J3n+fu06qQ+wRgmLu/4+6r3H0U4UO9y1peM8Tdv3P3KcBk4AV3n+Xu3xKOfrZLtqvQz1qGLkAT4Cp3X+7uLwNjgMNKbPMvd3/X3VcSCsC2Ff2hJb+obVViO9DdXyrjuTklbm8EvL/6jrt/b2bfEL7Bzy6x3fpAI2BiqAUAGOFbOoQjjWeqH5u2QH8zO63EY/WSnGX5ssTtZaXc3zC5XdGftTQbAXPcvbjEY58mr13tixK3lxIKhhQgFQDJZSWnqv2c8KELgJk1BpoD89Z4zdeED9NOHvoV1jQHKKvNuzJT484BLnf3yyvxmoqq6M9a1ms3NrNaJYrAb4GPajylpJ6agCQt7geOMbNtzaw+oVnnHXefXXKj5EPvTuB6M9sAwMxam9k+ySZ3Jfvpbma1kue2TJ77ktAmXxF3AieZ2U4WNDaz/c1sner9mED5P+vacr4D/ACca2Z1zawr0At4sAZySZ5RAZBUcPexwIXAY8B8wrf4Q8vY/C+EjtC3zew74CVgi2Q/7wLHANcTOoPH8fO37RuBP5rZIjO7qZw8Ewj9ADcDi5L3O7qKP96a+y7vZ70YGJWM8um7xmuXA72B/QhHQ7cCR1Wxn0PynLlrQRhJNzOrBawC2rr7Z7HziKSFjgAkH2wN/MgvOzdFpBwqAJJqZvZ/wCvAX5LmDxGpoIw2AZlZU2A44RuaA8cC04GHgHaEIW193X1RxkKIiEipMn0EcCPwnLtvCWwDTAUGA2PdvQPhTMjBGc4gIiKlyNgRgJmtC0wCNvESb5JMW9vV3eebWSvgVddMhiIiWZfJE8E2ARYAd5vZNsBE4AygpbvPB0iKwAbl7ahFixberl27KoWYsXAGHZp1qNJrY0lb5rTlhfRlTlteSF/mtOWF8jNPnDjxa3dfv6znM1kA6hAmszrN3d8xsxupRHOPmQ0ABgA0aN4guVV5K79eCS2q9tpY0pY5bXkhfZnTlhfSlzlteaECmU/k07XuIFOzzBHmNZld4v5uwNOETuBWyWOtgOnl7auoqMirqmhY1V8bS9oypy2ve/oypy2ve/oypy2ve/mZgQkeYzZQd/8CmGNmq9v3uxOmvB1NmDqX5PrJTGUQEZGyZXoyuNOA+5KFLmYRTsGvBTxsZscBnwF/ynAGEREpRUYLgLt/COxQylPdM/m+IiJSPp0JLCJSoFQAREQKlAqAiEiBUgEQESlQKgAiIgVKawKLZENxMSxYAJ9/Dl98AUuXwo8/hkudOtCoUbg0awZt2sCGG0LdurFTS55TARCpaYsWwRtvwAcfwH/+A5Mnw8yZsHJlxfdhFgpB587QuTM9FiyEXvOhVavM5ZaCowIgUl3Ll8Orr8KYMeF68mRwDx/im24KW28NffqED/TWrcO3+yZNoEEDqF8/FIZly+CHH+Cbb2Du3HCZOTMUkOef58qVK2H4RrDFFrDnntCrF3TvHl4vUkUqACJVsWIFPPcc3H8/PPMMfPcdNGwIu+4KffvC7rvDDjuEZp3qWr6cIy/ehnuaHw+vvAL33Qe33w7rrgsHHABHHgk9ekDt2tV/LykoKgAilTFtGgwfDvfeC19+CS1awB//CAceCHvtFYpATatXj6ntGsOAs+Hss+Gnn2DsWHjsMXjiiVCE2raF448Plw03rPkMkpc0CkikPO7hA7dnT+jYEW68EXbeGZ58MnTq3nVXaJLJxId/aerXD1nuuiu8/0MPwWabwYUXQrt2cNJJ8PHH2ckiqaYCIFIWd3j8cdhuu/DtfuJEuOQSmDcvPN67d/yROvXrhyanl16Cjz6Co4+Gu++GzTcPTUOzZ8fNJzlNBUBkTe7w7LOw445w8MGhg3bECPj00/Ate4NyF7GLo0OH0Dcwe3ZoKnr00dBpPGhQ6FwWWYMKgEhJH3wAe+wRmli++QZGjoQpU+CYY8KonTRo1QquvhpmzIB+/UKT1eabhyaj4uLY6SSHqACIQPiwP+WUMHJn2jS49VaYPh369w8naqVRmzbhQ//DD6FTp9BB3LUr/Pe/sZNJjlABkMLmHkb1bL453HEHDBwY2tJPPhnq1Yudrmb87nfh/IS77gpHM9tuC1ddBatWxU4mkakASOH69FPYe2844YTwIfnBB6G5pGnT2MlqXq1acOyx4ejmwAPhvPOgW7fwO5CCpQIghcedg19bEM7QffttuO22cILV734XO1nmrb9+GDY6alQoeJ07h/tSkFQApLAsWAD7789f7/sMfv/7MNXCSSeFaRsKhRkcdRRMmhSK4KGHhpFCK1bETiZZpgIgheOVV2CbbeDll7n60I3D2Pl27WKniqd9+/A7Oe00uP76MJ3El1/GTiVZpAIg+W/VKrj44jB52rrrwjvv8PCeGxTWt/6y1KsHN90E99wD774bRkFNmhQ7lWSJCoDkt4ULYb/94O9/D80eEyaEowD5pX79YPz4MCpqt93ghRdiJ5IsUAGQ/DVlSmjnHzcuDPUcOTJMwyyl23bb0Cnevn04EW7EiNiJJMNUACQ/PfEEdOkS5th/9VU47rjYidKhTRt4/fXQXHbccXD55bETSQapAEh+cYchQ+Cgg8LMnRMmhJk7peLWXTcsbtOvH1xwAfz1r+H3Knknpee4i5Ri1So4/fQwjcMhh4Qmn7TM35Nr6tYN5wo0agRXXhmOpK6/PpxQJnlDBUDyw9KlcPjhYY7+P/85THWgD6vqqVUrzC7aqBHccENYiOa22zR6Ko+oAEj6ffNNWBrxnXfCkMbTToudKH+YwdCh4UjqqqvCojdDh6oI5AkVAEm3L74IJzDNmBHmvz/44NiJ8o8ZXHFFOMq64YYwkurSS2OnkhqgAiDpNWdOWKlr7lx4+ukwckUywyx8+C9dCpddBo0bw+DBsVNJNakASDrNmhVms1y0KJy09Ic/xE6U/8xCn8DSpWE20WbNYMCA2KmkGlQAJH2mTQvf9n/8MSzWvsMOsRMVjtq1w+igxYvDmgmtW8P++8dOJVWkYRKSLjNnhm/+K1eGE7z04Z99deqEKaS32y4sSD9hQuxEUkUqAJIes2eHD/8VK+Dllwtj/v5c1aRJOFls/fXDCKxPPomdSKogowXAzGab2X/M7EMzm5A81szMXjSzGcn1epnMIHli3rzQ7LNkCbz4YljjVuLacEN49llYvjzMHfTtt7ETSSVl4whgT3ff1t1XH6sPBsa6ewdgbHJfpGxffhk+/BcsgOefD5OWSW7o2BH+9a/QNHfkkVBcHDuRVEKMJqA+wKjk9ijgwAgZJC0WLQpDPefMgWeeCbN7Sm7p2jUMEX3qKbjoothppBIyXQAceMHMJprZ6vFiLd19PkByvUGGM0haLVsGffrARx/B6NGw666xE0lZTjklzB562WXhhDxJBfMMzvJnZhu5++dmtgHwInAaMNrdm5bYZpG7/6ofICkYAwAaNG9Q1OmKqrX5Tv16Kh1bdKzSa2NJW+ZM5K1V7Fw1bBZdJy3mr8e356UdmtXo/vU7rnl1VxRz+9CP2HzuMo4evCVj6s/O+cwlpeF3vKbyMk88ceLEEs3vv+buWbkAFwPnANOBVsljrYDp5b22qKjIq6poWNVfG0vaMtd43uJi9xNPdAf3m26q2X0nCv53nCmff+7eqpX75pv7bjduGztNpaTmd1xCeZmBCb6Wz9aMNQGZWWMzW2f1bWBvYDIwGuifbNYfeDJTGSSlLrsMhg0LUw1oYrd0adUKHnwQZs7kr/d+qnUEclwm+wBaAm+Y2STgXeBpd38OuAroYWYzgB7JfZHg7rvhb38L6/decUXsNFIVu+8Ol17Kvu8tgjvvjJ1G1iJjU0G4+yzgV6tvu/s3gGbtkl8bNy7MLbPXXmENX005nF6DB/PWfUPY+fTTw8gtDd3NSToTWHLDzJlhKufNNoNHHgkrUkl61arFhce2g+bNw3QR338fO5GUQgVA4lu8GHr1CrefegqaNo2ZRmrI4nXqwv33h+I+aFDsOFIKFQCJa8WK8A3x44/DGaWbbRY7kdSkPfYIS3TeeWc4l0NyigqAxHXmmWFun2HDwoeF5J9LLoFttoHjjw/TekjOUAGQeO66C269Fc45B445JnYayZT69eG+++C770IR0NDQnKECIHG8916YPmCvveDKK2OnkUzr1CksKj9mjIaG5hAVAMm+r74KI35atYIHHggLjEj+O/30MKvr2WfDp5/GTiOoAEi2rVwJhx4KX38dOn1btIidSLKlVq1wfoc7nHiimoJygAqAZNd558Err4TFxbffPnYaybZ27UKT3/PPwz33xE5T8FQAJHseeQSuvRZOPRX69y9/e8lPp54Ku+wSRoB98UXsNAVNBUCy4+OPw3zxXbrA0KGx00hMtWqFEWA//KDJ/iJTAZDM++mncLJXnTphpsh69WInkti23DKsHvboo6EvSKJQAZDMO/dceP/9MNNn27ax00iu+POfwwlip58OS5bETlOQVAAksx5/HG66KbT39ukTO43kkrp14bbbYN48+PvfY6cpSCoAkjmzZ8Oxx8IOO8CQIbHTSC7aeedwdvANN8DkybHTFBwVAMmM5cvDeP/iYnjoIbX7S9muuirMAHvKKTo3IMtUACQz/vY3eOedMNpjk01ip5Fc1rx5KAKvv65zA7JMBUBq3rhxcPXVcMIJ8Mc/xk4jaXDssWGI8DnnwKJFsdMUDBUAqVmLF4f1fDfdVOP9peJq1Qodwt98AxdfHDtNwVABkJo1cGAY1XHvvdCkSew0kibbbhuOGm+9FaZNi52mIKgASM158MEw7/vf/gY77RQ7jaTRJZdAo0ahKUgyTgVAasacOXDyyaEd969/jZ1G0mqDDeCCC+Dpp+GFF2KnyXsqAFJtVuxhcrcVK0LTj+b3l+o4/fQwcmzQoDB9uGSMCoBU22EvfxWmeL7pptD5K1Id9evDNdfAlClaPSzDVACkeqZP55Qn5kGvXlrXV2rOQQfBHnuE/qTFi2OnyVsqAFJ1q1bBscfyU91aMGwYmMVOJPnCDK6/PgwL1ZrRGaMCIFV3000wfjzXHrJxWN9XpCZttx0ccUT4fzZvXuw0eUkFQKpmxoww2qdXL57dqVnsNJKvLrkkHGlecknsJHlJBUAqb9Wq0N7foEFY21dNP5Ip7dvDSSeFOaWmT4+dJu+oAEjl/eMf8Oab4dB8o41ip5F8d/754cvGhRfGTpJ3VACkclY3/RxwAPTrFzuNFIKWLeHss+GRR2DChNhp8ooKgFRccXGYq6V+fY36kew6+2xo0QLOOy92kryiAiAVd/fdYarna69V049k17rrhqagl16CsWNjp8kbKgBSMV9+GSbo2n33MHe7SLadfDK0aQMXXaSVw2pIxguAmdU2sw/MbExyv5mZvWhmM5Lr9TKdQWrAmWfC0qVq+pF46tcP/U9vvqmjgBqSjSOAM4CpJe4PBsa6ewdgbHJfctmzz4apns8/H7bcMnYaKWTHHgsbb6yjgBqS0QJgZm2A/YHhJR7uA4xKbo8CDsxkBqmmH34Ih94dO8Jf/hI7jRS61UcB48eH/gCplkwfAdwAnAsUl3ispbvPB0iuN8hwBqmOiy6CTz+FO+4If3wisR1zTDgKuPhiHQVUk3mGfoFmdgDQ091PMbOuwDnufoCZLXb3piW2W+Tuv+oHMLMBwACABs0bFHW6olOVckz9eiodW3Ss0mtjyZXMW3y2lH9eMZUnd23BFf3alrldruStjLRlTlteyGzm/xu3gPPu/4xTz+jAO1utWyP7zMff8cQTJ0509x3K3MDdM3IBrgTmArOBL4ClwL3AdKBVsk0rYHp5+yoqKvKqKhpW9dfGkhOZV650Lypy33BD90WL1rppTuStpLRlTlte9wxn/vFH9403dt9lF/fi4hrZZT7+joEJvpbP1ow1Abn7ee7ext3bAYcCL7t7P2A00D/ZrD/wZKYySDXceSdMnAhDh0LTprHTiPyS+gJqRIzzAK4CepjZDKBHcl9yyYIF4Y9rzz3h0ENjpxEp3THHQOvWcMUVsZOkVlYKgLu/6u4HJLe/cffu7t4huV6YjQxSCeedB0uWwM03a8y/5K769cMUEa++Cm+/HTtNKulMYPmlt98OU++edRZstVXsNCJrd8IJ0KyZVg2rIhUA+dmqVXDKKWGeH029K2nQpAmcdhqMHg2TJ8dOkzoqAPKzYcPggw9Cx+8668ROI1Ixp50GjRvDkCGxk6SOCoAEX30Vpnro1g369o2dRqTimjeHAQPggQfgk09ip0kVFQAJBg+G779Xx6+k06BBUKtWmKpcKkwFQOCtt8Jc/4MGhTl/RNKmTRs46qgwgOHLL2OnSQ0VgEJXXBzaUNXxK2l37rmwfDncckvsJKmhAlDo7rknnPE7ZEgYUSGSVptvDr16wa23wrJlsdOkggpAIfv++3DS1+9/D4cfHjuNSPUNGgTffBO+2Ei5VAAK2ZVXwvz5cOONoQNNJO123x223x6uvz40b8pa6a++UM2eDdddB0ccAV26xE4jUjPMwlHAtGnw3HOx0+Q8FYBCde654Vu/TqGXfNO3b5gkbujQ2ElyngpAIXrtNXjkkbDE48Ybx04jUrPq1oXTTw8Lx0+aFDtNTlMBKDSrVsGZZ4Zx03/+c+w0Iplxwglheojrr4+dJKepABSaUaPCfD9DhkCjRrHTiGTGeuvBscfC/feHgQ5SKhWAQrJkSVjopUsXOOyw2GlEMuuMM2DlSrjttthJcpYKQCG5+upwmvwNN2i+H8l/m24K++8Pd9wRzhCWX1EBKBSffx6GfR5yCOy0U+w0ItkxcGD40vPYY7GT5CQVgEJx0UXhcFjrp0oh6dEDNtsszHIrv6ICUAimTIERI+DUU2GTTWKnEcmeWrXC//vx48PgB/kFFYBCMHhwWOHrggtiJxHJvqOPDiPeNEvor6gA5LtXX4UxY8Kkb82bx04jkn1Nm0K/fnDffbBwYew0OUUFIJ8VF4eTvdq0CWdGihSqU0+FH38MCx/J/6gA5LNHHoEJE+Cyy6Bhw9hpROLp3Bl22y2sFaBZQv9HBSBf/fRTaPbp3Dkc/ooUuoEDYdYsePbZ2ElyhgpAvrr9dvjkk3DyV+3asdOIxHfQQbDhhuFvQwAVgPz07bdw6aWw116w996x04jkhrp1w/xAzzwDc+fGTpMTVADy0XXXhWXxhgzRlA8iJR1/fOgDGDEidpKcoAKQb776KiyE8ac/haXxRORn7duHo+Lhw8PU6AVOBSDfXHUVLFsGl1wSO4lIbhowAObMgeefj50kOhWAfDJnThjm1r8/bLll7DQiual3b2jZMswSWuBUAPLJpZeCe5j4TURKV7cuHHNMOEN+3rzYaaIqtwCY2Q5mdpaZXWNml5hZXzNrlo1wUgkzZoSOrZNOgrZtY6cRyW3HHx/6AAr8zOAyC4CZHW1m7wPnAQ2B6cBXwK7Ai2Y2ysx+m52YUq6LLoL69cOKXyKydptuGoZJF3hncJ21PNcY+IO7LyvtSTPbFugAfFbG8w2A14D6yfs86u4XJUcPDwHtgNlAX3dfVMX8AjBpEjzwQPjwb9kydhqRdBgwAPr2hRdfhH33jZ0mijKPANz9lrI+/JPnP3T3sWvZ909AN3ffBtgW2NfMugCDgbHu3gEYm9yX6rjwwjDj4TnnxE4ikh59+sD668Odd8ZOEs3ajgAAMLP2wGmEb+z/297de6/tde7uwPfJ3brJxYE+QNfk8VHAq8BfKpVafvbWW/DUU2Glr/XWi51GJD3q1QvzZN18czhxsgCVWwCAJ4C7gKeASk2jZ2a1gYnAZsAt7v6OmbV09/kA7j7fzDaoXGT5H/efm3003bNI5R19NFx/Pdx/f2isLjAWvqivZQOzd9y9WquIm1lT4HHCkcQb7t60xHOL3P1XX13NbAAwAKBB8wZFna7oVKX3nvr1VDq26Fil18ZS0cw7TPuO26+fwTWHbMxD3eLV0Xz+HeeKtOWF9GS+5/KpAOxwIqnIW1J5v+OJJ06c6O47lLmBu6/1AhwOXATsDGy/+lLe60rZz0XAOYTRRK2Sx1oB08t7bVFRkVdV0bCqvzaWCmUuLnbfdVf31q3dly3LfKi1yNvfcQ5JW173FGW+6SZ38EMu7Bg7SaWV9zsGJvhaPlsrciLY74ATgKuA65LLteW9yMzWT775Y2YNgb2AacBooH+yWX/gyQpkkDWNHQtvvBGagBo0iJ1GJL0OPxzq1qXX+MLrB6hIH8BBwCbuvryS+24FjEr6AWoBD7v7GDN7C3jYzI4jDCH9UyX3K+5w8cVhqcfjjoudRiTdmjeH3r3Z74UnYcWKcKZwgahIAZgENCWcBFZh7v5vYLtSHv8G6F6ZfckaXnoJ3nwzzPtTvwB7rkRq2tFH0+yxx8JqYb3XOsAxr1SkCaglMM3Mnjez0asvmQ4mZVg918/GG4fFLUSk+vbdl6/XrVNwU0NU5AhAM4vlkhdeCGP/b79d3/5FakqdOjyzU3OOGjMGFiwIJ4gVgLXNBWQA7j6utEvJbSRLVrf9//a3YTZDEakxY3ZpDitXwn33xY6SNWtrAnrFzE5bc8I3M6tnZt3MbBQ/j+aRbHj+eXj7bTj//HAWo4jUmFkbNYQdd4SRI2NHyZq1FYB9gVXAA2b2uZn918xmATOAw4Dr3X1kFjIK/Nz237ZtOHtRRGreUUeFyRUnT46dJCvWNhncj+5+q7v/AWhLGLmzvbu3dfcT3P3DbIUUwuiEd9+FCy7Qt3+RTOnbF2rXLphmoAqtCObuK9x9vrsvznAeKc3qtv927cJyjyKSGRtsAPvsE+YGKq7U1GeppCUh0+C55+C998K3/wI6SUUkiiOOgM8+C2fa5zkVgFznDpddFkb+HHVU7DQi+a9PH2jcuCCagSqyJvBAM9NE87GMGwfjx8Nf/qJv/yLZ0LgxHHggPPww/PRT7DQZVZEjgA2B98zsYTPbV2P/s+yyy2DDDXXWr0g29esHixeHwRd5rNwC4O4XENb+vQs4GphhZleY2aYZziZvvRVm/TznHM34KZJNe+0VOoTzvBmooqOAHPgiuawE1gMeNbOrM5hNLr88zFR44omxk4gUljp14JBDwnKr334bO03GVKQP4HQzmwhcDbwJ/M7dTwaKgP/LcL7C9cEH8PTTcNZZ0KRJ7DQihadfv9AH8NhjsZNkTEWOAFoAB7v7Pu7+iLuvAHD3YuCAjKYrZFdcAb/5DQwcGDuJSGHacUfYbLO8bgaqSB/A39z90zKem1rzkaT958vCt46BA0MREJHsMwvnBLzyCsybFztNRug8gBx0zHNfQMOGcOaZsaOIFLZDDw3n4jz6aOwkGaECkGs+/pi931sIJ58MLVrETiNS2LbcEjp3DucE5CEVgFwzZAirahmcfXbsJCICYTTQ+PEwZ07sJDVOBSCXzJkDI0fy5B9aQKtWsdOICIQZQgEeeSRujgxQAcgl11wD7vxzn5axk4jIapttBttvDw89FDtJjVMByBULFsDw4dCvH18011q/Ijmlb9+wHsfs2bGT1CgVgFxx882wbBmce27sJCKyptXNQHnWGawCkAt++CEUgD59oGPH2GlEZE3t24cTw1QApMYNHw4LF4Ypn0UkN/XtCxMnwscfx05SY1QAYluxAoYOhd12g513jp1GRMrypz+F6zw6ClABiO3BB8Pyc/r2L5Lb2raFLl1UAKSGuMPVV8PWW0PPnrHTiEh5DjkEPvwQPvoodpIaoQIQ0zPPwOTJYeSPFloTyX1//GO4zpNzAlQAYhoyJCz2fuihsZOISEW0aQO77AL/+lfsJDVCBSCWt96C11+HQYO02LtImhx8cGgGmjUrdpJqUwGIZcgQaNYMjj8+dhIRqYyDDgrXjz8eN0cNUAGIYepUePLJsOBL48ax04hIZWyyCWy7bV40A6kAxHDNNWHBl9NOi51ERKri4IPDFNHz58dOUi0ZKwBmtrGZvWJmU81sipmdkTzezMxeNLMZyfV6mcqQk+bOhXvvheOO04IvIml18MHh+oknosaorkweAawEznb3jkAX4FQz2woYDIx19w7A2OR+4bjxRigu1oIvImm21Vaw+eapbwbKWAFw9/nu/n5yewkwFWgN9AFGJZuNAg7MVIac8913cMcd4ZTydu1ipxGRqjILRwGvvBLm8UqprPQBmFk7YDvgHaClu8+HUCSADbKRISfcdVcoAoMGxU4iItV18MGwahU89VTsJFVm7p7ZNzBrAowDLnf3f5nZYndvWuL5Re7+q34AMxsADABo0LxBUacrOlXp/ad+PZWOLeJPsVx7lfP4hZP5Yr16DPjzFmvdNlcyV1Ta8kL6MqctL6Qvc6XzujPmvP8w/beNOPuUzTIXbC3KyzzxxIkT3X2HMjdw94xdgLrA88CgEo9NB1olt1sB08vbT1FRkVdV0bCqv7ZGPfSQO7g/8US5m+ZM5gpKW1739GVOW1739GWuUt7TT3evX999yZKaD1QB5WUGJvhaPlszOQrIgLuAqe4+tMRTo4H+ye3+wJOZypAz3OG668LaogccEDuNiNSUgw+Gn36CZ5+NnaRKMtkH8AfgSKCbmX2YXHoCVwE9zGwG0CO5n9/Gjw/riZ51FtSuHTuNiNSUXXeF9ddP7WigOpnasbu/AZQ1xWX3TL1vTrruujDtQ//+5W8rIulRu3ZYyvWhh8KRQP36sRNVis4EzrSZM8PJIiedpGkfRPJRnz6wZAmMGxc7SaWpAGTajTeG2T4HDoydREQyoXv3MLVLCoeDqgBk0sKFMGIEHH44tGoVO42IZELDhrD33jB6dBjwkSIqAJl0xx2wdGno/BWR/NWrV1jb+9//jp2kUlQAMmX5cvjHP6BHD+jcOXYaEcmkAw4I00OMHh07SaWoAGTKgw/C559r0jeRQtCyJey0kwqAENoBhw6FTp1C26CI5L/evWHChPDFLyVUADLh5Zdh0qQw6ZuVdSqEiOSV3r3D9ZgxcXNUggpAJtxwA2ywQRj9IyKFYautwnKRKWoGUgGoaTNnwtNPhxO/GjSInUZEssUsHAW89BL88EPsNBWiAlDTbr4Z6tQJBUBECkvv3mFKiBdfjJ2kQlQAatJ334UTv/r21YlfIoVo113hN79JTTOQCkBNGjkyzAlyxhmxk4hIDHXrQs+eoSN41arYacqlAlBTiovDiV9dusCOO8ZOIyKx9O4NCxaEKeBznApATXn22dABrG//IoVtn33CNNHPPBM7SblUAGrKjTdC69bwf/8XO4mIxLTeerDLLioABeO//w29/qecEtoARaSw7bcfvP8+zJ8fO8laqQDUhH/8I6wEdMIJsZOISC7o2TNcP/dc3BzlUAGorkWL4J//hCOOCGuDioh07gwbbZTzzUAqANU1fHiY8//002MnEZFcYRaOAl54AVasiJ2mTCoA1bFyZTjzd489YJttYqcRkVzSs2c4OXT8+NhJyqQCUB2jR4dVgDT0U0TW1L17GBSSw81AKgDVceON0K7dz9PAioistu66sNtuKgB56cMP4bXXYODAcNKHiMiaevaEyZNhzpzYSUqlAlBVN98MjRrBscfGTiIiuWq//cL1s8/GzVEGFYCqWLgQ7rsP+vULZ/2JiJSmY0do2zZnm4FUAKri7rvhxx/h1FNjJxGRXLZ6OOhLL4V1AnKMCkBlFRfDbbeFeb87d46dRkRyXc+eYYWw11+PneRXVAAq6/nn4eOP9e1fRCpmzz3DVDE52AykAlBZt9wCLVvCwQfHTiIiadC4cThZNAc7glUAKuOTT0IVHzAA6tWLnUZE0mKffWDatJwbDqoCUBm33Qa1asGJJ8ZOIiJpsvfe4fqFF+LmWIMKQEUtWwZ33QUHHhgWfhERqahOncLsoCoAKfXgg2H8/8CBsZOISNqYhaOAl17KqcXiM1YAzGyEmX1lZpNLPNbMzF40sxnJdTrOonIPnb+dOoXOHBGRytp77/Al8v33Yyf5n0weAYwE9l3jscHAWHfvAIxN7ue+d9+FiRPDko9msdOISBrttVe4fv75uDlKyFgBcPfXgIVrPNwHGJXcHgUcmKn3r1G33ALrrANHHhk7iYik1frrw/bb51Q/QLb7AFq6+3yA5HqDLL9/5S1YAA89BEcdFYqAiEhV7bMPvPVWWCgmB5i7Z27nZu2AMe6+dXJ/sbs3LfH8IncvtR/AzAYAAwAaNG9Q1OmKTlXKMPXrqXRs0bFKrwXo/9wXnPb4PP500VZ8slHDKu+nMqqbOdvSlhfSlzlteSF9mbORt2j6EoYN/YhBp2zKa9s0rfb+yss88cSJE919hzI3cPeMXYB2wOQS96cDrZLbrYDpFdlPUVGRV1XRsKq/1leudG/b1n3PPau+jyqoVuYI0pbXPX2Z05bXPX2Zs5L3xx/dGzd2P+WUGtldeZmBCb6Wz9ZsNwGNBvont/sDT2b5/Svn6afh008174+I1Iz69aFr15zpB8jkMNAHgLeALcxsrpkdB1wF9DCzGUCP5H7uuvXWcNJXnz6xk4hIvth7b5g5E2bNip2EOpnasbsfVsZT3TP1njVq1qwwXOvii6FOxn5NIlJo9tknXL/4YvRpZXQmcFnuuCOs9Xv88bGTiEg+2Xxz+O1vc6IZSAWgNMuXw4gR0KuX5v0RkZq1elqIsWNh5cqoUVQASvP442H8/0knxU4iIvlo773h22/hvfeixlABKM3tt0P79tCjR+wkIpKP9twzXI8dGzWGCsCapk2DV18Ni77U0q9HRDKgRQvYdlsVgJxzxx1Qty4cc0zsJCKSz7p3h/Hjw1ojkagAlLRsGYwcGdb7bdkydhoRyWfduoUBJ2++GS2CCkBJjz4KixZFH5srIgVg993DOUYRm4FUAEq6/fYwRrdr19hJRCTfNWkCO+0EL78cLYIKwGr//ndojzvpJC36IiLZ0b07TJgAixdHeXsVgNWGDQsTNfXvX/62IiI1oVs3KC6GceOivL0KAMD338M990DfvtCsWew0IlIounSBhg2jNQOpAAA8+CAsWaIzf0Uku+rXh912i9YRrAIAofP3d7+DnXeOnURECk23bjBlCnzxRdbfWgVgwgSYODEM/VTnr4hkW/dkhvwIzUAqAMOGQaNG0K9f7CQiUoi22w6aNlUByLpvv4X774fDD4ff/CZ2GhEpRLVrh8nhIvQDFHYBuPdeWLpUnb8iEle3bjB7dtaXiSzcAuAOd94J228PRUWx04hIIYvUD1C4BWDiRJg0CU44IXYSESl0W24JrVplvRmocAvA8OGh8/ewstauFxHJErPQDPTKK6F1IksKswB8/33o/O3bV52/IpIb9tgDvvwSPvooa29ZmAXgkUfCmb/HHx87iYhIsMce4frVV7P2loVZAO68Ezp2hF12iZ1ERCTo0CH0A2RxYrjCKwBTpsBbb4Vv/zrzV0RyhVk4Chg3Lmv9AIVXAIYPD2v+Hnlk7CQiIr/UtSt8/jnMnJmVtyusAvDTT/DPf8JBB8H668dOIyLyS6v7AbLUDFRYBeDxx2HhQnX+ikhu2mILaNkyax3BhVUAhg+Hdu1+PutORCSXZLkfoHAKwMcfh7PsjjsOahXOjy0iKbPHHjB3blbmBSqcT8IRI8IH/zHHxE4iIlK2rl3DdRb6AQqjAKxcCXffDT17QuvWsdOIiJStY8cwSEUFoIY88wzMn6/OXxHJfav7AbLQEVwYBeDOO8MZdvvvHzuJiEj59tgDPvssrBGQQXlfANZftDwcARx9NNSpEzuOiEj5sjQvUJQCYGb7mtl0M5tpZoMz+V69x38DxcVh9I+ISBp06gTNm2e8HyDrBcDMagO3APsBWwGHmdlWGXmz4mJ6v/l1mGd7000z8hYiIjWuVi3Yffe8PAL4PTDT3We5+3LgQaBPRt5p7Fhaf7Ncq36JSPp07Rr6AD77LGNvEaNRvDUwp8T9ucBOa25kZgOAAQANmjdghzt2qPQbXTLiE3ZqCL2+uZrld1xbxbjZN/XrqVX6eWNJW15IX+a05YX0Zc61vB0WLOWeWnDO1d14o3PTUrepdmZ3z+oF+BMwvMT9I4F/rO01RUVFXiVLlvhRg7es2msjKhpWxZ83krTldU9f5rTldU9f5pzLu2qV+7ffrnWT8jIDE3wtn60xmoDmAhuXuN8G+Dwj79SkCVPaN87IrkVEMqpWLVh33cy+RUb3Xrr3gA5m1t7M6gGHAqMj5BARKWhZ7wNw95VmNhB4HqgNjHD3KdnOISJS6KKcGeXuzwDPxHhvEREJ8v5MYBERKZ0KgIhIgVIBEBEpUCoAIiIFSgVARKRAqQCIiBQoFQARkQJlYbqI3GZmC4BPq/jyFsDXNRgnG9KWOW15IX2Z05YX0pc5bXmh/Mxt3X39sp5MRQGoDjOb4O65M8VfBaQtc9ryQvoypy0vpC9z2vJC9TOrCUhEpECpAIiIFKhCKAB3xA5QBWnLnLa8kL7MacsL6cuctrxQzcx53wcgIiKlK4QjABERKUVeFwAz29fMppvZTDMbHDvPmsxsYzN7xcymmtkUMzsjebyZmb1oZjOS6/ViZy3JzGqb2QdmNia5n+t5m5rZo2Y2Lfld75yCzGcl/ycmm9kDZtYglzKb2Qgz+8rMJpd4rMx8ZnZe8nc43cz2yaHM1yT/L/5tZo+bWdNcyVxa3hLPnWNmbmYtSjxW6bx5WwDMrDZwC7AfsBVwmJltFTfVr6wEznb3jkAX4NQk42BgrLt3AMYm93PJGcDUEvdzPe+NwHPuviWwDSF7zmY2s9bA6cAO7r41YeGkQ8mtzCOBfdd4rNR8yf/pQ4FOyWtuTf4+s20kv878IrC1u3cGPgLOg5zJPJJf58XMNgZ6AJ+VeKxKefO2AAC/B2a6+yx3Xw48CPSJnOkX3H2+u7+f3F5C+GBqTcg5KtlsFHBglIClMLM2wP7A8BIP53LedYHdgbsA3H25uy8mhzMn6gANzawO0IiwbnbOZHb314CFazxcVr4+wIPu/pO7fwLMJPx9ZlVpmd39BXdfmdx9m7BGOeRA5jJ+xwDXA+cCJTtwq5Q3nwtAa2BOiftzk8dykpm1A7YD3gFauvt8CEUC2CBitDXdQPjPV1zisVzOuwmwALg7abYabmaNyeHM7j4PuJbwDW8+8K27v0AOZ06UlS8tf4vHAs8mt3Mys5n1Bua5+6Q1nqpS3nwuAFbKYzk55MnMmgCPAWe6+3ex85TFzA4AvnL3ibGzVEIdYHvgNnffDviBHGruKU3Sdt4HaA9sBDQ2s35xU1VLzv8tmtn5hCbZ+1Y/VMpmUTObWSPgfOBvpT1dymPl5s3nAjAX2LjE/TaEw+icYmZ1CR/+97n7v5KHvzSzVsnzrYCvYuVbwx+A3mY2m9Ck1s3M7iV380L4fzDX3d9J7j9KKAi5nHkv4BN3X+DuK4B/AbuQ25mh7Hw5/bdoZv2BA4Aj/Odx8bmYeVPCl4JJyd9gG+B9M9uQKubN5wLwHtDBzNqbWT1CB8noyJl+wcyM0DY91d2HlnhqNNA/ud0feDLb2Urj7ue5ext3b0f4fb7s7v3I0bwA7v4FMMfMtkge6g78lxzOTGj66WJmjZL/I90J/UO5nBnKzjcaONTM6ptZe6AD8G6EfL9iZvsCfwF6u/vSEk/lXGZ3/4+7b+Du7ZK/wbnA9sn/8arldfe8vQA9CT37HwPnx85TSr5dCYdp/wY+TC49geaEURQzkutmsbOWkr0rMCa5ndN5gW2BCcnv+QlgvRRk/jswDZgM3APUz6XMwAOE/okVyQfRcWvLR2i6+BiYDuyXQ5lnEtrOV//93Z4rmUvLu8bzs4EW1cmrM4FFRApUPjcBiYjIWqgAiIgUKBUAEZECpQIgIlKgVABERAqUCoCISIFSARCpJDNraGbjKjM7pJkNNLNjMplLpLJ0HoBIJZnZqUAdd7+xEq9pBLzpYT4ikZygIwCRhJntmCwM0sDMGicLsmxdyqZHkExzYGZdk6OBh83sIzO7ysyOMLN3zew/ZrYpgIdpBmabWdanQRYpS53YAURyhbu/Z2ajgcuAhsC97v6L1ZiSeaU2cffZJR7eBuhImLt9FjDc3X9vYYW304Azk+0mALuRI/PgiKgAiPzSJYSJBH8krMq1phbA4jUee8+TefDN7GPgheTx/wB7ltjuK2DLmgwrUh1qAhL5pWZAE2AdoEEpzy8r5fGfStwuLnG/mF9+yWqQvF4kJ6gAiPzSHcCFhIVBhqz5pLsvAmqbWWnFoTybE2b3FMkJKgAiCTM7Cljp7vcDVwE7mlm3UjZ9gTCVd2X9AXipGhFFapSGgYpUkpltBwxy9yMz+RqRTNMRgEglufsHwCuVORGM0Hl8YYYiiVSJjgBERAqUjgBERAqUCoCISIFSARARKVAqACIiBUoFQESkQP0/FNx9IVwwOyMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "#\n",
    "#debug = True\n",
    "debug = False\n",
    "#\n",
    "def drag(cd, area, rho, velx, vely):\n",
    "    '''\n",
    "    Return horizontal and vertical drag force on body given its drag coefficient, area, \n",
    "    density of medium in which it's moving and horizontal and vertical velocity.\n",
    "    '''\n",
    "    v2 = velx**2 + vely**2\n",
    "    sinTheta = vely/np.sqrt(v2)\n",
    "    cosTheta = velx/np.sqrt(v2)\n",
    "    Dx = -0.5*cd*rho*area*v2*cosTheta\n",
    "    Dy = -0.5*cd*rho*area*v2*sinTheta\n",
    "    return Dx, Dy\n",
    "#\n",
    "# Initialise variables. Approximate (maximum) number of steps in calculation\n",
    "if debug:\n",
    "    nSteps = 5\n",
    "else:\n",
    "    nSteps = 1000\n",
    "#\n",
    "# Arrays for x and y coordinates\n",
    "xStep, yStep = np.zeros(nSteps), np.zeros(nSteps) # m\n",
    "#\n",
    "# Initial values of x and y coordinates\n",
    "xE, yE = 0.0, 0.0  # m\n",
    "#\n",
    "# Initial values of velocities in x and y dorections\n",
    "ux, uy = 30, 40 # m/s\n",
    "vX, vY = ux, uy # m/s\n",
    "#\n",
    "# Acceleration due to gravity\n",
    "g = -9.81 # m/s**2\n",
    "#\n",
    "# Time for projectile fight without air resistance\n",
    "tMax = -2*uy/g # s\n",
    "#\n",
    "# Time interval for each step\n",
    "dt = tMax/nSteps\n",
    "#\n",
    "# Drag coefficient\n",
    "CD = 0.47 \n",
    "#\n",
    "# Radius of spherical projectile\n",
    "rad = 0.025 # m\n",
    "#\n",
    "# Frontal area of projectile\n",
    "Area = np.pi*rad**2 # m**2\n",
    "#\n",
    "# Densities of air and of projectile\n",
    "rhoAir = 1.2 # kg/m**3\n",
    "rhoProj = 2000.0 # kg/m**3\n",
    "#\n",
    "# Mass of projectile\n",
    "mProj = 4/3*np.pi*rad**3*rhoProj # kg\n",
    "print(\" \")\n",
    "print(\"Radius {:5.3f} m, mass {:5.3f} kg.\".format(rad, mProj))\n",
    "print(\"Initial position ({:5.2f}, {:5.2f}) m.\".format(xE, yE))\n",
    "print(\"Initial velocity ({:5.2f}, {:5.2f}) m/s.\".format(vX, vY))\n",
    "print(\"Maximum time {:5.3f} s, time step {:5.3e} s.\".format(tMax, dt))\n",
    "#\n",
    "# Step number\n",
    "iStep = 0\n",
    "#\n",
    "# Calculate while projectile above ground level (or initial step) and number of steps less than maximum\n",
    "# Look at values used to control while loop\n",
    "if debug:\n",
    "    print(\" \")\n",
    "    print(\"Initial tests:\")\n",
    "    print(f\"ye = {yE} and iStep = {iStep}\")\n",
    "    print(f\"yE > 0  = {yE > 0 }, iStep == 0 = {iStep == 0} and iStep < nSteps = {iStep < nSteps}.\")\n",
    "    print(f\"Hence yE > 0 or iStep == 0 = {yE > 0 or iStep == 0}\")\n",
    "    print(f\"And (yE > 0 or iStep == 0) and iStep < nSteps = {(yE > 0 or iStep == 0) and iStep < nSteps}\")\n",
    "#\n",
    "while (yE > 0 or iStep == 0) and iStep < nSteps:\n",
    "    #\n",
    "    # Save coordinate values\n",
    "    xStep[iStep], yStep[iStep] = xE, yE\n",
    "    #\n",
    "    # Update coordinates for this step\n",
    "    xE = xE + vX*dt\n",
    "    yE = yE + vY*dt\n",
    "    #\n",
    "    # Find components of drag force\n",
    "    dragX, dragY = drag(CD, Area, rhoAir, vX, vY)\n",
    "    #\n",
    "    # Update velocities in x and y directions\n",
    "    vX = vX + dragX/mProj*dt\n",
    "    vY = vY + dragY/mProj*dt + g*dt\n",
    "    #\n",
    "    # Increment the step number\n",
    "    iStep = iStep + 1\n",
    "    if debug:\n",
    "        print(\" \")\n",
    "        print(\"Tests for loop {iStep}\")\n",
    "        print(f\"ye = {yE} and iStep = {iStep}\")\n",
    "        print(f\"yE > 0  = {yE > 0 }, iStep == 0 = {iStep == 0} and iStep < nSteps = {iStep < nSteps}.\")\n",
    "        print(f\"Hence yE > 0 or iStep == 0 = {yE > 0 or iStep == 0}\")\n",
    "        print(f\"And (yE > 0 or iStep == 0) and iStep < nSteps = {(yE > 0 or iStep == 0) and iStep < nSteps}\")\n",
    "#\n",
    "print(\" \")\n",
    "print(\"Final step is number {:d}, actual flight time {:5.3f} s.\".format(iStep, iStep*dt))\n",
    "print(\" \")\n",
    "plt.figure(figsize = (6, 5))\n",
    "plt.title(\"Projectile motion\")\n",
    "plt.xlabel(\"x (m)\")\n",
    "plt.ylabel(\"y (m)\")\n",
    "plt.plot(xStep[0:iStep], yStep[0:iStep], linestyle = '-', color = 'r')\n",
    "plt.grid(color = 'g')\n",
    "plt.show()"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}