Lecture.04-Projectile-Motion.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Projectile Motion: The Trajectory of a Cannon Shell"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Going forward from the bicyle problem, we turn our attention to two-dimensional projectile motion.  Our specific example will be the shell shot by a large cannon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "\n",
    "YouTubeVideo(\"hlW6hZkgmkA\", playlist=\"hlW6hZkgmkA\", autoplay=1, loop=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In the absence of air-resistance, the equations of motion (Newton's 2nd Law) are:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\frac{d^2 x}{dt^2} =& 0 \\\\\n",
    "\\frac{d^2 y}{dt^2} =& -g\n",
    "\\end{eqnarray}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "These are second-order differential equations.  A standard technique for solving higher order differential equations is to rewrite them as a system of first order differential equations.  That is,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\frac{d x}{dt} =& v_x \\\\\n",
    "\\frac{d v_x}{dt} =& 0 \\\\\n",
    "\\frac{d y}{dt} =& v_y \\\\\n",
    "\\frac{d v_y}{dt} =& -g\n",
    "\\end{eqnarray}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "There are now four equations written out for the variables $x, y, v_x, v_y$.  As before, we can solve this problem by using Euler's method:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "x_{i+1} &= x_i + v_{x,i} \\Delta t \\\\\n",
    "v_{x,i+1} &= v_{x,i} \\\\\n",
    "y_{i+1} &= y_i + v_{y,i} \\Delta t \\\\\n",
    "v_{y,i+1} &= v_{y,i} - g \\Delta t\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This algorithm allows us to estimate the position and velocity of a projectile given the intial values of $x, y, v_x$, and $v_y$.\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Like in the bicycle problem, air resistance is important. To introduce air resistance we can model the drag as\n",
    "\n",
    "$$ F_{drag} = - B_2 v^2 $$\n",
    "\n",
    "like we did in the bicycle problem where $v$ is the magnitude of the velocity vector. Since force is vector so we need to split this into horizontal and vertical components."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This force is a vector and always directed in the opposite direction of the velocity vector.\n",
    "\n",
    "\\begin{align}\n",
    "F_{drag,x} = F_{drag} \\cos \\theta = F_{drag} (v_x / v) \\\\\n",
    "F_{drag,y} = F_{drag} \\sin \\theta = F_{drag} (v_y / v)\n",
    "\\end{align}\n",
    "\n",
    "we can write\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "F_{drag, x} = - B_2 v v_x \\\\\n",
    "F_{drag, y} = - B_2 v v_y \n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We can add these drag forces to our algorithm like we did last lecture."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "\\begin{align}\n",
    "x_{i+1} &= x_i + v_{x,i} \\Delta t \\\\\n",
    "v_{x,i+1} &= v_{x,i} - \\frac{B_2 v v_x}{m} \\Delta t \\\\\n",
    "y_{i+1} &= y_i + v_{y,i} \\Delta t \\\\\n",
    "v_{y,i+1} &= v_{y,i} - g \\Delta t - \\frac{B_2 v v_x}{m} \\Delta t\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Simulation of position vs. time for a cannon ball"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Initialize variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# acceleration due to gravity, m/s^2\n",
    "g = 9.81\n",
    "# time step, s\n",
    "dt = 0.1\n",
    "# max time, s\n",
    "tmax = 200"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In this problem, the variables are `x`, `y`, `vx`, and `vy`.  In the code below we create numpy arrays for `x` and `y` since we want to record the values of position for over many different timesteps.  \n",
    "\n",
    "The velocity components `vx` and `vy` are simple variables that represent only the value of the velocity at the *current* time step.\n",
    "\n",
    "Finally, in this solution we don't even have time, `t`, listed explicited.  There is still a time step, `dt`, and a maximum time, `tmax`, but no array for time.\n",
    "\n",
    "Depending on the problem we are solving, we can choose what variables we want to record. If you preferred, you can have arrays for `vx`, `vy`, and `t` as well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Compute steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def update(i, x, y, vx, vy, Bm):\n",
    "    # update positions\n",
    "    x[i+1] = x[i] + vx*dt\n",
    "    y[i+1] = y[i] + vy*dt\n",
    "        \n",
    "    # calculate velocity magnitude\n",
    "    v = np.sqrt(vx**2 + vy**2)\n",
    "        \n",
    "    # compute air drag\n",
    "    Fdrag_x = -(Bm)*v*vx\n",
    "    Fdrag_y = -(Bm)*v*vy\n",
    "        \n",
    "    # update velocities\n",
    "    vx_new = vx + Fdrag_x * dt\n",
    "    vy_new = vy + (-g + Fdrag_y) * dt\n",
    "    \n",
    "    return vx_new, vy_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def solve(theta=45, v0=700, Bm=4e-5):\n",
    "    \"\"\"\n",
    "    Given \n",
    "        theta, firing angle (degrees)\n",
    "        v0, initial velocity (m/s)\n",
    "        Bm, drag coefficient per unit mass (1/m)\n",
    "    Compute\n",
    "        x,y trajectory of a cannon shell\n",
    "    \"\"\"\n",
    "    \n",
    "    # Allocate memory for arrays and set to zero\n",
    "    N = round(tmax / dt)\n",
    "    x = np.zeros(N)\n",
    "    y = np.zeros(N)\n",
    "    \n",
    "    # Initial values\n",
    "    x[0] = 0\n",
    "    y[0] = 0\n",
    "    vx = v0*np.cos(np.deg2rad(theta))\n",
    "    vy = v0*np.sin(np.deg2rad(theta))\n",
    "\n",
    "    # Calculate the solution\n",
    "    for i in range(N-1):\n",
    "        vx, vy = update(i, x, y, vx, vy, Bm)\n",
    "        \n",
    "        # stop when y <= 0\n",
    "        if y[i+1] <= 0:\n",
    "            break\n",
    "\n",
    "    # truncate arrays to when shell has hit the ground (y<=0)\n",
    "    N = i+2\n",
    "    x = x[:N]\n",
    "    y = y[:N]\n",
    "\n",
    "    # interpolate position if shell is below ground level\n",
    "    if y[i+1] < 0:\n",
    "        r = -y[i]/y[i+1]\n",
    "        x[i+1] = (x[i]+r*x[i+1])/(r+1)\n",
    "        y[i+1] = 0 \n",
    "\n",
    "    return x, y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the main `for` loop has an `if` statement inside. This code checks if the projectile has hit the ground.  The truncation and interpolation steps make it so the trajectory ends exactly at the ground. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Plotting routine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def plot(x, y, label='', style='-'):\n",
    "    plt.title('Trajectory of cannon shell')\n",
    "    # show results in km\n",
    "    plt.plot(x/1000, y/1000, style, label=label)\n",
    "    plt.xlabel('x (km)')\n",
    "    plt.ylabel('y (km)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice we are plotting `x` vs `y`.  For 2D projectile problems this is the best way of visualing the solution. \n",
    "\n",
    "If you wanted, you could also plot `x` vs `t` and `y` vs `t` as two separate plots.  In that case, you would need to calculate and store the time `t` as an array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Driver code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots()\n",
    "\n",
    "x, y = solve(theta=theta, Bm=4e-5)\n",
    "plot(x, y, label='With air drag')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Exercise\n",
    "\n",
    "> 1. Run the again model for without air drag. How does this change the range?\n",
    "> 2. Run the model for several different firing angles. For what angle is the range maximized?\n",
    "> 3. Compare your angle of maximum range with and without air drag. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Effect of Air Density\n",
    "\n",
    "Clearly air drag plays an imporant role in this problem.  However, the plots showed that the shells travel many kilometers up into the air.  The air density will be lower up high in the atmosphere than down at sea level. \n",
    "\n",
    "From last lecture we know that air drag is proportional to the density of air. How important will this high altitude effect be on the on the trajectory of the cannon shell?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The details of how air density varies as a function of altitude will be covered in Phys 3400: Thermodynamics.  Here, let's assume that the atmosphere is an *isothermal* (constant temperature) ideal gas.  The pressure $p$ depends on alitude like\n",
    "\n",
    "$$ p(y) = p(0) \\exp(-mgy/k_B T) $$\n",
    "\n",
    "where $m$ is the average mass of an air molecule, $y$ is the height from some reference point like sea level, $k_B$ is the Boltzmann's constant, and $T$ is the absolute temperature.  For an ideal gas, the pressure so this leads to\n",
    "\n",
    "$$ \\rho = \\rho_0 \\exp(-y/y_0) $$\n",
    "\n",
    "where $y_0 = k_B T / m g \\approx 1.0 \\times 10^4$ m, and $\\rho_0$ is the density at sea level ($y=0$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The air drag force, which is proportional to density, then can be modified to be a function of altitude\n",
    "\n",
    "$$ F_{drag}^* = \\frac{\\rho}{\\rho_0} F_{drag}(y=0) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This value $y_0$ is sometimes called the scale height.  Let's modify our computation steps to include the effect of this scale height."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def update(i, x, y, vx, vy, Bm, y0=None):\n",
    "    \n",
    "    # update positions\n",
    "    x[i+1] = x[i] + vx*dt\n",
    "    y[i+1] = y[i] + vy*dt\n",
    "        \n",
    "    # calculate velocity magnitude\n",
    "    v = np.sqrt(vx**2 + vy**2)\n",
    "        \n",
    "    # compute air drag\n",
    "    Fdrag_x = -(Bm)*v*vx\n",
    "    Fdrag_y = -(Bm)*v*vy\n",
    "    \n",
    "    # correct for density change\n",
    "    if y0 is not None:\n",
    "        rho_rho0 = np.exp(-y[i]/y0)\n",
    "        Fdrag_x = rho_rho0 * Fdrag_x\n",
    "        Fdrag_y = rho_rho0 * Fdrag_y\n",
    "    \n",
    "    # update velocities\n",
    "    vx = vx + Fdrag_x * dt\n",
    "    vy = vy + (-g + Fdrag_y) * dt\n",
    "    \n",
    "    return vx, vy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def solve(theta=45, v0=700, Bm=4e-5, y0=None):\n",
    "    \"\"\"\n",
    "    Given \n",
    "        theta, firing angle (degrees)\n",
    "        v0, initial velocity (m/s)\n",
    "        Bm, drag coefficient per unit mass (1/m)\n",
    "        y0, scale height of atmosphere\n",
    "    Compute\n",
    "        x,y trajectory of a cannon shell\n",
    "    \"\"\"\n",
    "    \n",
    "    # Allocate memory for arrays and set to zero\n",
    "    N = round(tmax / dt)\n",
    "    x = np.zeros(N)\n",
    "    y = np.zeros(N)\n",
    "    \n",
    "    # Initial values\n",
    "    x[0] = 0\n",
    "    y[0] = 0\n",
    "    vx = v0*np.cos(np.deg2rad(theta))\n",
    "    vy = v0*np.sin(np.deg2rad(theta))\n",
    "\n",
    "    # Calculate the solution\n",
    "    for i in range(N-1):\n",
    "        vx, vy = update(i, x, y, vx, vy, Bm, y0)\n",
    "        \n",
    "        # stop when y <= 0\n",
    "        if y[i+1] <= 0:\n",
    "            break\n",
    "\n",
    "    # truncate arrays to when shell has hit the ground (y<=0)\n",
    "    N = i+2\n",
    "    x = x[:N]\n",
    "    y = y[:N]\n",
    "\n",
    "    # interpolate position if shell is below ground level\n",
    "    if y[i+1] < 0:\n",
    "        r = -y[i]/y[i+1]\n",
    "        x[i+1]=(x[i]+r*x[i+1])/(r+1)\n",
    "        y[i+1]=0 \n",
    "\n",
    "    return x, y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots()\n",
    "\n",
    "x, y = solve(theta=45, Bm=4e-5, y0=None)\n",
    "plot(x, y, style=':', label='Without density correction')\n",
    "x, y = solve(theta=35, Bm=4e-5, y0=None)\n",
    "plot(x, y, style=':', label='Without density correction')\n",
    "x, y = solve(theta=45, Bm=4e-5, y0=1e4)\n",
    "plot(x, y, label='With density correction')\n",
    "x, y = solve(theta=35, Bm=4e-5, y0=1e4)\n",
    "plot(x, y, label='With density correction')\n",
    "plt.ylim(0,12)\n",
    "plt.legend(loc='upper left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Exercise\n",
    "> In our model of the cannon shell trajectory we have assumed that the acceleration due to gravity, $g$, is constant. It will also depend on altitude. Add this to the model and calculate how much it affects the range."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The following formula approximates the Earth's gravity variation with altitude:\n",
    "\n",
    "$$\n",
    "g_{h}=g_{0}\\left({\\frac  {r_{{\\mathrm  {e}}}}{r_{{\\mathrm  {e}}}+h}}\\right)^{2}$$\n",
    "\n",
    "where\n",
    "\n",
    " * $g_h$ is the gravitational acceleration at height h above sea level.\n",
    " * $r_e$ is the [Earth's mean radius](https://en.wikipedia.org/wiki/Earth_radius)\n",
    " * $g_0$ is the [standard gravitational acceleration](https://en.wikipedia.org/wiki/Standard_gravitational_acceleration)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}