{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Case study: Throwing Axe\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Throwing axe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our favorite event at Lumberjack Competitions is axe throwing.  The axes used for this event typically weigh 1.5 to 2 kg, with handles roughly 0.7 m long.  They are thrown overhead at a target typically 6 m away and 1.5 m off the ground.  Normally, the axe makes one full rotation in the air to hit the target blade first, with the handle close to vertical.\n",
    "\n",
    "![Diagram of throwing axe](diagrams/throwingaxe1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a version of `make_system` that sets the initial conditions.\n",
    "\n",
    "The state variables are x, y, theta, vx, vy, omega, where theta is the orientation (angle) of the axe in radians and omega is the angular velocity in radians per second.\n",
    "\n",
    "I chose initial conditions based on videos of axe throwing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = UNITS.meter\n",
    "s = UNITS.second\n",
    "kg = UNITS.kilogram\n",
    "radian = UNITS.radian\n",
    "\n",
    "def make_system():\n",
    "    \"\"\"Makes a System object for the given conditions.\n",
    "    \n",
    "    returns: System with init, ...\n",
    "    \"\"\"\n",
    "    P = Vector(0, 2) * m\n",
    "    V = Vector(8, 4) * m/s\n",
    "    \n",
    "    init = State(x=P.x, y=P.y, theta=2, \n",
    "                 vx=V.x, vy=V.y, omega=-7)\n",
    "\n",
    "    t_end = 1.0 * s\n",
    "    \n",
    "    return System(init=init, t_end=t_end,\n",
    "                  g = 9.8 * m/s**2,\n",
    "                  mass = 1.5 * kg,\n",
    "                  length = 0.7 * m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make a `System`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = make_system()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "system.init"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a simple starting place, I ignore drag, so `vx` and `omega` are constant, and `ay` is just `-g`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def slope_func(state, t, system):\n",
    "    \"\"\"Computes derivatives of the state variables.\n",
    "    \n",
    "    state: State (x, y, x velocity, y velocity)\n",
    "    t: time\n",
    "    system: System object with length0, m, k\n",
    "    \n",
    "    returns: sequence (vx, vy, ax, ay)\n",
    "    \"\"\"\n",
    "    x, y, theta, vx, vy, omega = state\n",
    "    unpack(system)\n",
    "\n",
    "    ax = 0\n",
    "    ay = -g\n",
    "    alpha = 0\n",
    "\n",
    "    return vx, vy, omega, ax, ay, alpha"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As always, let's test the slope function with the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope_func(system.init, 0, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "results, details = run_ode_solver(system, slope_func, max_step=0.05)\n",
    "details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing the results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simplest way to visualize the results is to plot the state variables as a function of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(results.x, label='x')\n",
    "plot(results.y, label='y')\n",
    "\n",
    "decorate(xlabel='Time (s)',\n",
    "         ylabel='Position (m)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(results.theta, label='theta', color='C2')\n",
    "\n",
    "decorate(xlabel='Time (s)',\n",
    "         ylabel='Angle (radian)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the velocities the same way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(results.vx, label='vx')\n",
    "plot(results.vy, label='vy')\n",
    "\n",
    "decorate(xlabel='Time (s)',\n",
    "         ylabel='Velocity (m/s)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(results.omega, label='omega', color='C2')\n",
    "\n",
    "decorate(xlabel='Time (s)',\n",
    "         ylabel='Angular velocity (rad/s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to visualize the results is to plot y versus x.  The result is the trajectory through the plane of motion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(results.x, results.y, label='trajectory', color='C3')\n",
    "\n",
    "decorate(xlabel='x position (m)',\n",
    "         ylabel='y position (m)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "**Exercise:**  Find the starting conditions that make the final height of the COG as close as possible to 1.5 m.  Ideally, the final angle should be a little past vertical."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  Compute the total velocity of the leading edge of the axe at the point of impact, that is, the sum of velocity due to translation and rotation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_frame(theta):\n",
    "    rhat = Vector(pol2cart(theta, 1))\n",
    "    that = rhat.perp()\n",
    "    return rhat, that"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# quick, get those state variables into Vector objects!\n",
    "state = get_last_value(results)\n",
    "x, y, theta, vx, vy, omega = state\n",
    "P = Vector(x, y)\n",
    "V = Vector(vx, vy)\n",
    "rhat, that = make_frame(theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Animation\n",
    "\n",
    "NOTE: This section needs to be updated.\n",
    "\n",
    "Animating this system is a little more complicated, if we want to show the shape and orientation of the axe.\n",
    "\n",
    "It is useful to construct a frame with $\\hat{r}$ along the handle of the axe and $\\hat{\\theta}$ perpendicular."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're ready to animate the results.  The following figure shows the frame and the labeled points A, B, C, and D.\n",
    "\n",
    "![Diagram of the axe with reference frame](diagrams/throwingaxe2.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "decorate(xlabel='x position (m)',\n",
    "         ylabel='y position (m)',\n",
    "         xlim=[0, 8.1],\n",
    "         ylim=[0, 5.5],\n",
    "         legend=False)\n",
    "\n",
    "\n",
    "\n",
    "for t, state in system.results.iterrows():\n",
    "    x, y, theta, vx, vy, omega = state\n",
    "    P = Vector(x, y)\n",
    "    rhat, that = make_frame(theta)\n",
    "    \n",
    "    # plot the handle\n",
    "    A = P - l1 * rhat\n",
    "    B = P + l2 * rhat\n",
    "    plot_segment(A, B, color='red', update=True)\n",
    "\n",
    "    # plot the axe head\n",
    "    C = B + l2 * that\n",
    "    D = B - l2 * that\n",
    "    plot_segment(C, D, color='black', linewidth=10, update=True)\n",
    "\n",
    "    # plot the COG\n",
    "    plot(x, y, 'bo', update=True)\n",
    "    sleep(0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "During the animation, the parts of the axe seem to slide around relative to each other.  I think that's because the lines and circles get rounded off to the nearest pixel.\n",
    "\n",
    "Here's the final state of the axe at the point of impact (assuming the target is 8 m away)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}