{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Clock Pendulum with Air Drag and Joint Friction\n",
    "\n",
    "This notebook builds on the previous one by introducing both a nonlinear pendulum and a nonlinear damping effect through [Coulomb friction](https://en.wikipedia.org/wiki/Friction#Dry_friction). Students will be able to work with both a linear and nonlinear version of the same system (a clock pendulum) in order to compare the free response in both cases.\n",
    "\n",
    "After the completion of this assignment students will be able to:\n",
    "\n",
    "- compute the free response of a non-linear compound pendulum\n",
    "- estimate the period of oscillation of a nonlinear system\n",
    "- compute the free response of a non-linear compound pendulum with Coulomb damping\n",
    "- compare the non-linear behavior to the linear behavior\n",
    "- identify the function that governs the decay envelope for Coulomb damping\n",
    "\n",
    "![](fig/05/coulomb-pendulum.png)\n",
    "\n",
    "A common source of damping and energy dissipation is friction, which is a specific type of damping. With viscous damping, we had a torque that was linearly proportional to the angular velocity:\n",
    "\n",
    "$$T_f = c l^2 \\dot{\\theta}$$\n",
    "\n",
    "This simple source of energy dissipation is a reasonable mathematical model for many phenomena, but it isn't often a good model for dry friction between to hard materials. One very useful model of friction is [Coulomb friction](https://en.wikipedia.org/wiki/Coulomb_damping). Coulomb friction behaves more like:\n",
    "\n",
    "$$T_f = \\begin{cases} -\\frac{2}{3}\\mu R F_N & \\dot{\\theta} > 0 \\\\ 0 & \\dot{\\theta} = 0 \\\\ \\frac{2}{3}\\mu R F_N & \\dot{\\theta} < 0 \\end{cases}$$\n",
    "\n",
    "where $\\mu$ is a coefficient of sliding friction, $R$ is the outer radius of the joint contact (assuming disc/disc contact, see http://www.iitg.ernet.in/kd/Lecture%20Notes/ME101-Lecture15-KD_DivI.pdf for an explanation), and $F_N$ is the normal force in the pivot joint. Here the damping torque is constant, always impeding the motion of the system.\n",
    "\n",
    "This can also be more simply written using the [signum function](https://en.wikipedia.org/wiki/Sign_function):\n",
    "\n",
    "$$ T_f = -\\frac{2}{3}\\mu R F_n \\text{sgn}\\left( \\dot{\\theta} \\right)$$\n",
    "\n",
    "To start import some of the common packages we will need:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear vs Nonlinear\n",
    "\n",
    "`resonance` has both a version of the clock pendulum with the viscous damping (a model of air drag) and Coulomb friction (a model of joint friction). The first is a [linear system](https://en.wikipedia.org/wiki/Linear_system) due to the torque being linear in the angular velocity and the second is a [nonlinear system](https://en.wikipedia.org/wiki/Nonlinear_system) because the torque is nonlinear with respect to the velocity. Import them both:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from resonance.nonlinear_systems import ClockPendulumSystem as NonLinPendulum\n",
    "from resonance.linear_systems import ClockPendulumSystem as LinPendulum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "non_lin_sys = NonLinPendulum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(NonLinPendulum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "non_lin_sys.constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lin_sys = LinPendulum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(LinPendulum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lin_sys.constants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "**Exercise**\n",
    "\n",
    "Simulate each system for 5 seconds with an initial angle of 1.0 degrees and plot the trajectory of the angle from each response on a single plot to see if there are any differences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "initial_angle = np.deg2rad(1.0)\n",
    "non_lin_sys.coordinates['angle'] = initial_angle\n",
    "lin_sys.coordinates['angle'] = initial_angle\n",
    "\n",
    "duration = 5.0\n",
    "non_lin_traj = non_lin_sys.free_response(duration)\n",
    "lin_traj = lin_sys.free_response(duration)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write your answer here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "\n",
    "ax.plot(non_lin_traj.index, non_lin_traj.angle)\n",
    "ax.plot(lin_traj.index, lin_traj.angle)\n",
    "\n",
    "ax.legend(['Nonlinear', 'Linear']);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "**Exercise**\n",
    "\n",
    "Create a function called `compare_lin_to_nonlin` that accepts a single angle value in degrees as the only argument and creates a plot just like the one you made above. Don't forget axis labels and a legend so that we can tell the two lines apart. Include the initial angle in the title of the plot. It should look something like:\n",
    "\n",
    "```python\n",
    "def compare_lin_to_nonlin(initial_angle):\n",
    "    # write your code here (hint copy most of it from above and modify)\n",
    "    fig, ax = plt.subplots(1, 1)\n",
    "    ax.plot(non_lin_traj.index, np.rad2deg(non_lin_traj.angle))\n",
    "    ax.plot(lin_traj.index, np.rad2deg(lin_traj.angle))\n",
    "    # write the legend, labels, title, etc here\n",
    "    ax.grid()\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "def compare_lin_to_nonlin(initial_angle):\n",
    "    initial_angle = np.deg2rad(initial_angle)\n",
    "    non_lin_sys.coordinates['angle'] = initial_angle\n",
    "    lin_sys.coordinates['angle'] = initial_angle\n",
    "    non_lin_traj = non_lin_sys.free_response(duration)\n",
    "    lin_traj = lin_sys.free_response(duration)\n",
    "    fig, ax = plt.subplots(1, 1)\n",
    "    ax.plot(non_lin_traj.index, np.rad2deg(non_lin_traj.angle))\n",
    "    ax.plot(lin_traj.index, np.rad2deg(lin_traj.angle))\n",
    "    ax.legend(['Nonlinear', 'Linear'])\n",
    "    ax.set_xlabel('Time [s]')\n",
    "    ax.set_ylabel('Angle [deg]')\n",
    "    ax.set_title('Initial angle: {:1.0f} deg'.format(np.rad2deg(initial_angle)))\n",
    "    ax.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "**Exercise**\n",
    "\n",
    "Try out some angles from 0 to 180 degrees and view the graphs to see if there is anything interesting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "compare_lin_to_nonlin(1)\n",
    "compare_lin_to_nonlin(25)\n",
    "compare_lin_to_nonlin(87)\n",
    "compare_lin_to_nonlin(130)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "**Exercise**\n",
    "\n",
    "Make a plot of initial angle versus period for the nonlinear pendulum. Report on your what you learn from this plot. The code should look something like:\n",
    "\n",
    "```python\n",
    "periods = []\n",
    "angles = np.deg2rad(np.linspace(0, 90))\n",
    "for angle in angles:\n",
    "    # fill in the loop here\n",
    "    \n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(np.rad2deg(angles), periods)\n",
    "ax.set_ylabel('Period [s]')\n",
    "ax.set_xlabel('Initial Angle [deg]')\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "from resonance.functions import estimate_period\n",
    "periods = []\n",
    "angles = np.deg2rad(np.linspace(0, 90))\n",
    "for angle in angles:\n",
    "    non_lin_sys.coordinates['angle'] = angle\n",
    "    traj = non_lin_sys.free_response(5.0)\n",
    "    periods.append(estimate_period(traj.index, traj.angle))\n",
    "    \n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(np.rad2deg(angles), periods)\n",
    "ax.set_ylabel('Period [s]')\n",
    "ax.set_xlabel('Initial Angle [deg]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adding Damping\n",
    "\n",
    "Below we add a small amount of viscous damping to the linear pendulum and a small amount of Coulomb friction to the nonlinear pendulum. If you compare the trajectories of the angle at an initial angle of 5 degrees you see some differences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lin_sys.constants['viscous_damping'] = 0.1  # Ns/m\n",
    "non_lin_sys.constants['coeff_of_friction'] = 0.1  # Nm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare_lin_to_nonlin(5.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Curve Fitting an Exponential Decay Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "non_lin_traj = non_lin_sys.free_response(5.0, sample_rate=500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "**Exercise**\n",
    "\n",
    "Use your curve fitting function for $\\theta(t) = A e^{\\lambda t} \\cos(\\omega t + \\phi)$ and see how well it fits the nonlinear trajectory.\n",
    "\n",
    "```python\n",
    "def exp_decayed_oscillation(times, amplitude, decay_constant, frequency, phase_shift):\n",
    "    # write your function here and return the answer\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "def exp_decayed_oscillation(times, amplitude, decay_constant, frequency, phase_shift):\n",
    "    return amplitude * np.exp(decay_constant * times) * np.cos(frequency * times + phase_shift)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since there is something funny going on at the end of the simulation, only fit to the first 3 seconds of data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "popt, pcov = curve_fit(exp_decayed_oscillation,\n",
    "                       non_lin_traj[:3.0].index, non_lin_traj[:3.0].angle,\n",
    "                       p0=(5.0, -0.001, 2 * np.pi, 0.0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(non_lin_traj.index, non_lin_traj.angle, '.')\n",
    "best_fit_angle = exp_decayed_oscillation(non_lin_traj.index, popt[0], popt[1], popt[2], popt[3])\n",
    "ax.plot(non_lin_traj.index, best_fit_angle)\n",
    "ax.legend(['Data', 'Fit'])\n",
    "ax.set_xlabel('Time [s]')\n",
    "ax.set_ylabel('Angle [rad]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "**Exercise**\n",
    "\n",
    "Now try a function that looks like:\n",
    "\n",
    "$$ \\theta(t) = (m t + A) \\cos(\\omega t + \\phi) $$\n",
    "\n",
    "This is a linear decaying function instead of an exponentially decaying function.\n",
    "\n",
    "```python\n",
    "def lin_decayed_oscillation(times, m, A, frequency, phi):\n",
    "    # write function here and return result\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "def lin_decayed_oscillation(times, m, A, frequency, phi):\n",
    "    return (m * times + A) * np.cos(frequency * times + phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write your answer here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "popt, pcov = curve_fit(lin_decayed_oscillation,\n",
    "                       non_lin_traj[:3.0].index, non_lin_traj[:3.0].angle,\n",
    "                       p0=(-0.1, 0.07, 2 * np.pi, 0.0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(non_lin_traj.index, non_lin_traj.angle, '.')\n",
    "best_fit_angle = lin_decayed_oscillation(non_lin_traj.index, popt[0], popt[1], popt[2], popt[3])\n",
    "ax.plot(non_lin_traj.index, best_fit_angle)\n",
    "ax.legend(['Data', 'Fit'])\n",
    "ax.set_xlabel('Time [s]')\n",
    "ax.set_ylabel('Angle [rad]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Investigate the stick-slip area\n",
    "\n",
    "Notice that there is a point in time where the pendulum stops oscillating. At this point the pendulum does not have enough kinetic energy to overcome the friction. We can look more closely at this by simulating at a higher sample rate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "**Exercise**\n",
    "\n",
    "Add a measurement to the system that computes the friction torque based on the defintion at the beginning of the notebook. Assume that $T_N=1$ Newton-meter, $\\mu=0.2$, and $\\theta_0=2.0$ deg. Simulate the system and plot the angle, angular velocity, and friction torque in 3 subplots that share the same X axis (time)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "def calculate_friction(coeff_of_friction, angle_vel):\n",
    "    return -coeff_of_friction * np.sign(angle_vel)\n",
    "non_lin_sys.add_measurement('friction', calculate_friction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "non_lin_sys.constants['coeff_of_friction'] = 0.2  # Nm\n",
    "non_lin_sys.coordinates['angle'] = np.deg2rad(2.0)\n",
    "non_lin_traj = non_lin_sys.free_response(2.0, sample_rate=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(3, 1, sharex=True)\n",
    "\n",
    "axes[0].plot(non_lin_traj.index,\n",
    "        np.rad2deg(non_lin_traj.angle))\n",
    "axes[0].set_ylabel('Angle [deg]')\n",
    "axes[0].grid()\n",
    "\n",
    "axes[1].plot(non_lin_traj.index,\n",
    "        np.rad2deg(non_lin_traj.angle_vel))\n",
    "axes[1].set_ylabel('Angular Velocity [deg/s]')\n",
    "axes[1].grid()\n",
    "\n",
    "axes[2].plot(non_lin_traj.index,\n",
    "             non_lin_traj.friction)\n",
    "axes[2].set_ylabel('Friction Torque [Nm]')\n",
    "axes[2].set_xlabel('Time [s]')\n",
    "axes[2].grid()\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "non_lin_traj.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.rad2deg(non_lin_traj.angle).tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualization\n",
    "\n",
    "Finally, visualize the animation with different friction coefficients to see the behavior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "non_lin_sys.coordinates['angle'] = np.deg2rad(5.0)\n",
    "non_lin_traj = non_lin_sys.free_response(4.0, sample_rate=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ani = non_lin_sys.animate_configuration(interval=1000/100)  # interval should be in milliseconds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import HTML, display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "html = ani.to_jshtml()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(HTML(html))"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}