diff --git a/Lab.02-Skydiver-Problem.ipynb b/Lab.02-Skydiver-Problem.ipynb new file mode 100644 index 0000000..5f22701 --- /dev/null +++ b/Lab.02-Skydiver-Problem.ipynb @@ -0,0 +1,437 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Skydiver Problem\n", + "## Lab 2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Name:\n", + "### Student ID:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Instructions:\n", + "*Complete the notebook by reading all of the markdown cells and filling in the *blanks*. In the code cells, you will need to replace with the correct code everywhere you see a `___`. You will also need to add new code as appropirate. Please complete this lab during the lab slot and ask for help from the instructor or the TA as often as needed.*" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This lab involves modelling a skydiver falling through the air. Initially, air resistance is very small, but upon deployment of a parachute, air drag suddenly increases to a large value. The mass of the jumper is 60 kg, acceleration due to gravity is 9.81 m/s^2 , the jumper starts from a height of 1000m and opens the parachute at about 200 m above the ground. The jumper feels a drag force $F=-k v |v|$, where $k$ changes to a high value after the parachute is opened." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Model the action of the parachute with a hyperbolic tangent" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " \n", + " To model the action of the parachute we will use the hyperbolic tangent function. Plot `np.tanh(x-5)` from x=0 to 10. Don't worry about making the graph look perfect, but make sure the range of the graph goes from -1.2 to 1.2." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x = np.arange(___, ___, 0.01)\n", + "y = ___\n", + "\n", + "fig, axes = plt.subplots()\n", + "plt.ylim(___, ___)\n", + "plt.plot(x, y)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As you can see from the graph, to a good approximation the function goes from one constant value to another. Now plot the function 6 - 2 tanh( (x-5)/0.1 ) from $x$ = 0 to 10. Choose an appropriate range so that you can see both the top and bottom of your plot." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# your code to make plot goes here\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice how we have transformed the function so that it changes more abruptly, has a wider range and different end values. Let's now choose appropriate parameters to model our system, and store them in variables\n", + "\n", + "* m = 60 - mass of jumper in kg\n", + "* g = 9.81 - acceleration due to gravity\n", + "* h = 1000 - initial height in m\n", + "* y0 = 200 - height at which the jumper opens the parachute\n", + "* kd = 0.48 - drag coefficient without parachute in kg /m\n", + "* kp = 30.24 - drag coefficient with parachute fully deployed\n", + "* d = 20 - roughly the distance (m) the jumper falls during the parachute opening" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# initialize all the variables here\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define a new Python function that implements the the formula\n", + "$$k(y) = kp - (kp-kd)/2 (1+tanh( 2 (y-y0)/d )).$$ \n", + "\n", + "k(y) is our mathematical model for the drag coefficient as a function of the jumper's height above the ground." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def k(y):\n", + " return ___" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the function k(y), with y ranging from 0 to 1000." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "y = np.arange(0, 1000, 0.1)\n", + "\n", + "fig, axes = plt.subplots()\n", + "plt.xlabel(\"height (m)\")\n", + "plt.ylabel(\"drag coefficient (kg/m)\")\n", + "plt.plot(y, k(y))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Check your values for k(0) and k(1000). They should be 30.24 and 0.48 respectively. \n", + "\n", + "k(0) \n", + "k(1000)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(k(0), k(1000))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Solving the problem numerically" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " Now that we've modelled the skydiver problem with a height-varying drag coefficient, let's find the height and velocity as functions of time using Euler's method, assuming initial velocity of zero. \n", + "\n", + "The equation that we want to solve is\n", + "$$ \\frac{d^2y}{dt^2} = -g - k(y)/m \\frac{dy}{dt} \\left| \\frac{dy}{dt} \\right| \\;, $$\n", + "\n", + "along with appropriate initial conditions.\n", + "\n", + "The first task is to rewrite this second-order differential equation as pair of two first-order differential equations by introducing $v = \\frac{dy}{dt}$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\\begin{align}\n", + "\\frac{dy}{dt} &= \\mbox{___} \\\\\n", + "\\frac{dv}{dt} &= \\mbox{___}\n", + "\\end{align}\n", + "\n", + "*Determine what should be on the right hand side and edit this Markdown cell to replace the `\\mbox{___}` accordingly.*" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Implement your numerical solution below using Euler's method." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Computation routine" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def compute():\n", + " tmax = 400\n", + " dt = ___ # you should check that your choice of dt gives a good solution\n", + " \n", + " N = round(tmax/dt)\n", + " \n", + " # allocate memory for y, v, and t arrays\n", + " y = np.zeros(N)\n", + " v = np.zeros(N)\n", + " t = np.zeros(N)\n", + " \n", + " # initialize variables\n", + " y[0] = ___\n", + " v[0] = ___\n", + " t[0] = ___\n", + " \n", + " # solve using Euler's method\n", + " for i in range(N-1):\n", + " y[i+1] = ___\n", + " v[i+1] = ___\n", + " t[i+1] = ___\n", + " \n", + " if y[i+1] <= 0:\n", + " break\n", + " \n", + " # truncate arrays to when skydiver has landed (y<=0)\n", + " N = i+2\n", + " y = y[:N]\n", + " v = v[:N]\n", + " t = t[:N]\n", + " \n", + " return (y, v, t)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Run your numerical model and plot the height $y$ vs time $t$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "y, v, t = compute()\n", + "\n", + "# plot the solution\n", + "fig, axes = plt.subplots()\n", + "\n", + "plt.plot(___, ___)\n", + "plt.xlabel(___)\n", + "plt.ylabel(___)\n", + "plt.title('Skydiver height vs time')\n", + "plt.ylim(0, h)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3. Analysis of solution" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Estimate the total time the skydiver spends in the air.\n", + "\n", + "* Easier (approximate) method: Use the plot to estimate time when when y = 0.\n", + "* Harder (exact) method: Use linear interpolation to estimate the time when time y = 0 (we did this is a previous lecture).\n", + "* *Try either approach (or both!) as you prefer.*" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "t_landing = ___\n", + "print(\"The skydiver lands after\", t_landing, \"s.\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the velocity $v$ versus time $t$. Remember to use good axes labels with units and add a title!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot the velocity vs time\n", + "fig, axes = plt.subplots()\n", + "\n", + "plt.plot(___, ___)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "What is the approximate terminal speed of the skydiver before the parachute is deployed? What about after?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# You can look at the actually values of v (using `print(v[i])` at some appropirately chosen values of `i`)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Answer: Terminal speed without the parachute is `___` m/s, while with the parachute the it is `___` m/s." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Calculate and plot the acceleration of the skydiver (include axes labels with units).\n", + "\n", + "* *Hint: you will need to calculate acceleration as the change in velocity divided by the change in time. You have already computed both velocity and time.*" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "a = np.zeros_like(v)\n", + "N = len(a)\n", + "for i in range(N-1):\n", + " a[i+1] = ___\n", + " \n", + "fig, axes = plt.subplots()\n", + "plt.plot(t,a)\n", + "plt.xlabel(___)\n", + "plt.ylabel(___)\n", + "plt.ylim(___, ___)\n", + "plt.title(___)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "What is the largest acceleration felt by the skydiver during his descent (express your answer as a multiple of g)?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# the max() function returns the maximum of an array\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Answer: The largest acceleration felt is approximately `___` g." + ] + } + ], + "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 +} diff --git a/Lecture.04-Projectile-Motion.ipynb b/Lecture.04-Projectile-Motion.ipynb new file mode 100644 index 0000000..abc4188 --- /dev/null +++ b/Lecture.04-Projectile-Motion.ipynb @@ -0,0 +1,710 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "< [Bicycle Racing: The Effect of Air Resistance](Lecture.03-Bicyle-Racing.ipynb) | [Contents](Index.ipynb) | [Baseball](Lecture.05-Baseball.ipynb) >" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "plt.style.use('default')" + ] + }, + { + "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)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "< [Bicycle Racing: The Effect of Air Resistance](Lecture.03-Bicyle-Racing.ipynb) | [Contents](Index.ipynb) | [Baseball](Lecture.05-Baseball.ipynb) >" + ] + } + ], + "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.3" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +}