{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bungee Dunk Revisited"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Modeling and Simulation in Python*\n",
    "\n",
    "Copyright 2021 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# install Pint if necessary\n",
    "\n",
    "try:\n",
    "    import pint\n",
    "except ImportError:\n",
    "    !pip install pint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# download modsim.py if necessary\n",
    "\n",
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "    \n",
    "download('https://github.com/AllenDowney/ModSimPy/raw/master/modsim.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import functions from modsim\n",
    "\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous case study, we simulated a bungee jump with a model that took into account gravity, air resistance, and the spring force of the bungee cord, but we ignored the weight of the cord.\n",
    "\n",
    "It is tempting to say that the weight of the cord doesn't matter because it falls along with the jumper.  But that intuition is incorrect, as explained by [Heck, Uylings, and Kędzierska](http://iopscience.iop.org/article/10.1088/0031-9120/45/1/007).  As the cord falls, it transfers energy to the jumper.   They derive a differential equation that relates the acceleration of the jumper to position and velocity:\n",
    "\n",
    "$a = g + \\frac{\\mu v^2/2}{\\mu(L+y) + 2L}$ \n",
    "\n",
    "where $a$ is the net acceleration of the jumper, $g$ is acceleration due to gravity, $v$ is the velocity of the jumper, $y$ is the position of the jumper relative to the starting point (usually negative), $L$ is the length of the cord, and $\\mu$ is the mass ratio of the cord and jumper.\n",
    "\n",
    "If you don't believe this model is correct, [this video might convince you](https://www.youtube.com/watch?v=X-QFAB0gEtE).\n",
    "\n",
    "Following the previous case study, we'll model the jump with the following assumptions:\n",
    "\n",
    "1. Initially the bungee cord hangs from a crane with the attachment point 80 m above a cup of tea.\n",
    "\n",
    "2. Until the cord is fully extended, it applies a force to the jumper as explained above.\n",
    "\n",
    "3. After the cord is fully extended, it obeys [Hooke's Law](https://en.wikipedia.org/wiki/Hooke%27s_law); that is, it applies a force to the jumper proportional to the extension of the cord beyond its resting length.\n",
    "\n",
    "4. The jumper is subject to drag force proportional to the square of their velocity, in the opposite of their direction of motion.\n",
    "\n",
    "First I'll create a `Param` object to contain the quantities we'll need:\n",
    "\n",
    "1. Let's assume that the jumper's mass is 75 kg and the cord's mass is also 75 kg, so `mu=1`.\n",
    "\n",
    "2. The jumpers's frontal area is 1 square meter, and terminal velocity is 60 m/s.  I'll use these values to back out the coefficient of drag.\n",
    "\n",
    "3. The length of the bungee cord is `L = 25 m`.\n",
    "\n",
    "4. The spring constant of the cord is `k = 40 N / m` when the cord is stretched, and 0 when it's compressed.\n",
    "\n",
    "I adopt the coordinate system and most of the variable names from [Heck, Uylings, and Kędzierska](http://iopscience.iop.org/article/10.1088/0031-9120/45/1/007).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = Params(y_attach = 80,   # m,\n",
    "                 v_init = 0,     # m / s,\n",
    "                 g = 9.8,        # m/s**2,\n",
    "                 M = 75,         # kg,\n",
    "                 m_cord = 75,    # kg\n",
    "                 area = 1,       # m**2,\n",
    "                 rho = 1.2,      # kg/m**3,\n",
    "                 v_term = 60,    # m / s,\n",
    "                 L = 25,         # m,\n",
    "                 k = 40,         # N / m\n",
    "               )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now here's a version of `make_system` that takes a `Params` object as a parameter.\n",
    "\n",
    "`make_system` uses the given value of `v_term` to compute the drag coefficient `C_d`.\n",
    "\n",
    "It also computes `mu` and the initial `State` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(params):\n",
    "    \"\"\"Makes a System object for the given params.\n",
    "    \n",
    "    params: Params object\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    M, m_cord = params.M, params.m_cord\n",
    "    g, rho, area =  params.g, params.rho, params.area\n",
    "    v_init, v_term = params.v_init, params.v_term\n",
    "    \n",
    "    # back out the coefficient of drag\n",
    "    C_d = 2 * M * g / (rho * area * v_term**2)\n",
    "    \n",
    "    mu = m_cord / M\n",
    "    init = State(y=params.y_attach, v=v_init)\n",
    "    t_end = 8\n",
    "\n",
    "    return System(params, C_d=C_d, mu=mu,\n",
    "                  init=init, t_end=t_end)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make a `System`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "system1 = make_system(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`drag_force` computes drag as a function of velocity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def drag_force(v, system):\n",
    "    \"\"\"Computes drag force in the opposite direction of `v`.\n",
    "    \n",
    "    v: velocity\n",
    "    \n",
    "    returns: drag force in N\n",
    "    \"\"\"\n",
    "    rho, C_d, area = system.rho, system.C_d, system.area\n",
    "\n",
    "    f_drag = -np.sign(v) * rho * v**2 * C_d * area / 2\n",
    "    return f_drag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's drag force at 20 m/s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "drag_force(20, system1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function computes the acceleration of the jumper due to tension in the cord.\n",
    "\n",
    "$a_{cord} = \\frac{\\mu v^2/2}{\\mu(L+y) + 2L}$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cord_acc(y, v, system):\n",
    "    \"\"\"Computes the force of the bungee cord on the jumper:\n",
    "    \n",
    "    y: height of the jumper\n",
    "    v: velocity of the jumpter\n",
    "    \n",
    "    returns: acceleration in m/s\n",
    "    \"\"\"\n",
    "    L, mu = system.L, system.mu\n",
    "    \n",
    "    a_cord = -v**2 / 2 / (2*L/mu + (L+y))\n",
    "    return a_cord"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's acceleration due to tension in the cord if we're going 20 m/s after falling 20 m."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = -20\n",
    "v = -20\n",
    "cord_acc(y, v, system1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now here's the slope function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_func1(t, state, system):\n",
    "    \"\"\"Compute derivatives of the state.\n",
    "    \n",
    "    state: position, velocity\n",
    "    t: time\n",
    "    system: System object containing g, rho,\n",
    "            C_d, area, and mass\n",
    "    \n",
    "    returns: derivatives of y and v\n",
    "    \"\"\"\n",
    "    y, v = state\n",
    "    M, g = system.M, system.g\n",
    "    \n",
    "    a_drag = drag_force(v, system) / M\n",
    "    a_cord = cord_acc(y, v, system)\n",
    "    dvdt = -g + a_cord + a_drag\n",
    "    \n",
    "    return v, dvdt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As always, let's test the slope function with the initial params."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope_func1(0, system1.init, system1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll need an event function to stop the simulation when we get to the end of the cord."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def event_func1(t, state, system):\n",
    "    \"\"\"Run until y=-L.\n",
    "    \n",
    "    state: position, velocity\n",
    "    t: time\n",
    "    system: System object containing g, rho,\n",
    "            C_d, area, and mass\n",
    "    \n",
    "    returns: difference between y and y_attach-L\n",
    "    \"\"\"\n",
    "    y, v = state   \n",
    "    return y - (system.y_attach - system.L)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can test it with the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "event_func1(0, system1.init, system1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "results1, details1 = run_solve_ivp(system1, slope_func1, \n",
    "                                  events=event_func1)\n",
    "details1.message"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the plot of position as a function of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_position(results, **options):\n",
    "    results.y.plot(**options)\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Position (m)')\n",
    "    \n",
    "plot_position(results1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use `min` to find the lowest point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "min(results1.y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, Phase 1 ends when the jumper reaches an altitude of 55 m."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Phase 2\n",
    "\n",
    "Once the jumper has falled more than the length of the cord, acceleration due to energy transfer from the cord stops abruptly.  As the cord stretches, it starts to exert a spring force.  So let's simulate this second phase."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`spring_force` computes the force of the cord on the jumper:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def spring_force(y, system):\n",
    "    \"\"\"Computes the force of the bungee cord on the jumper:\n",
    "    \n",
    "    y: height of the jumper\n",
    "    \n",
    "    Uses these variables from system:\n",
    "    y_attach: height of the attachment point\n",
    "    L: resting length of the cord\n",
    "    k: spring constant of the cord\n",
    "    \n",
    "    returns: force in N\n",
    "    \"\"\"\n",
    "    L, k = system.L, system.k\n",
    "    \n",
    "    distance_fallen = system.y_attach - y\n",
    "    extension = distance_fallen - L\n",
    "    f_spring = k * extension\n",
    "    return f_spring"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spring force is 0 until the cord is fully extended.  When it is extended 1 m, the spring force is 40 N. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "spring_force(55, system1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "spring_force(56, system1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The slope function for Phase 2 includes the spring force, and drops the acceleration due to the cord."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_func2(t, state, system):\n",
    "    \"\"\"Compute derivatives of the state.\n",
    "    \n",
    "    state: position, velocity\n",
    "    t: time\n",
    "    system: System object containing g, rho,\n",
    "            C_d, area, and mass\n",
    "    \n",
    "    returns: derivatives of y and v\n",
    "    \"\"\"\n",
    "    y, v = state\n",
    "    M, g = system.M, system.g\n",
    "    \n",
    "    a_drag = drag_force(v, system) / M\n",
    "    a_spring = spring_force(y, system) / M\n",
    "    dvdt = -g + a_drag + a_spring\n",
    "    \n",
    "    return v, dvdt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial state for Phase 2 is the final state from Phase 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_final = results1.index[-1]\n",
    "t_final"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "state_final = results1.iloc[-1]\n",
    "state_final"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And that gives me the starting conditions for Phase 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "system2 = System(system1, t_0=t_final, init=state_final)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we run Phase 2, setting the direction of the event function so it doesn't stop the simulation immediately. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "results2, details2 = run_solve_ivp(system2, slope_func2)\n",
    "details2.message"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_final = results2.index[-1]\n",
    "t_final"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the results on the same axes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_position(results1, label='Phase 1')\n",
    "plot_position(results2, label='Phase 2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And get the lowest position from Phase 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "min(results2.y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see how big the effect of the cord is, I'll collect the previous code in a function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_two_phases(params):\n",
    "    system1 = make_system(params)\n",
    "    results1, details1 = run_solve_ivp(system1, slope_func1, \n",
    "                                       events=event_func1)\n",
    "    t_final = results1.index[-1]\n",
    "    state_final = results1.iloc[-1]\n",
    "    \n",
    "    system2 = system1.set(t_0=t_final, init=state_final)\n",
    "    results2, details2 = run_solve_ivp(system2, slope_func2)\n",
    "    return pd.concat([results1, results2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run both phases and get the results in a single `TimeFrame`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = run_two_phases(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_position(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "params_no_cord = params.set(m_cord=1)\n",
    "results_no_cord = run_two_phases(params_no_cord);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plot_position(results, label='m_cord = 75 kg')\n",
    "plot_position(results_no_cord, label='m_cord = 1 kg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "min(results_no_cord.y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = min(results.y) - min(results_no_cord.y)\n",
    "diff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The difference is about a meter, which could certainly be the difference between a successful bungee dunk and a bad day."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}