{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. Vertical Vibration of Quarter Car Model\n",
    "\n",
    "This notebook introduces the base excitation system by examning the behavior of a quarter car model.\n",
    "\n",
    "After the completion of this assignment students will be able to:\n",
    "\n",
    "- excite a system with a sinusoidal input\n",
    "- understand the difference in transient and steady state solutions\n",
    "- create a frequency response plot\n",
    "- define resonance and determine the parameters that cause resonance\n",
    "\n",
    "![](fig/07/quarter-car.jpg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib notebook\n",
    "from resonance.linear_systems import SimpleQuarterCarSystem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sys = SimpleQuarterCarSystem()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simple quarter car model has a suspension stiffness and damping, along with the sprung car mass in kilograms, and a travel speed parameter in meters per second."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.speeds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "# A sinusoidal road\n",
    "\n",
    "The road is described as:\n",
    "\n",
    "$$y(t) = Ysin\\omega_b t$$\n",
    "\n",
    "where $Y$ is the amplitude of the sinusoidal road undulations and $\\omega_b$ is the frequency of the a function of the car's speed. If the distance between the peaks (amplitude 0.01 meters) of the sinusoidal road is 6 meters and the car is traveling at 7.5 m/s calculate what the frequency will be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "Y = 0.01  # m\n",
    "v = sys.constants['travel_speed']\n",
    "bump_distance = 6  # m\n",
    "wb = v / bump_distance * 2 * np.pi  # rad /s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(wb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "Now with the amplitude and frequency set you can use the `sinusoidal_base_displacing_response()` function to simulate the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "traj = sys.sinusoidal_base_displacing_response(Y, wb, 20.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj.plot(subplots=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've written an animation for you. You can play it with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.animate_configuration(fps=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "Try different travel speeds and see what kind of behavior you can observe. Make sure to set the `travel_speed` constant and the frequency value for `sinusoidal_base_displacing_response()` to be consistent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "# Transmissibility\n",
    "\n",
    "When designing a car the designer wants the riders to feel comfortable and to isolate them from the road's bumps. There are two important aspects to investigate. The first is called *displacement transmissibility* and is a ratio between the ampitude of the steady state motion and the ampitude of the sinusoidal base displacement. So in our case this would be:\n",
    "\n",
    "$$ \\frac{X}{Y}(\\omega_b) = \\frac{\\textrm{Steady State Amplitude}}{\\textrm{Base Displacement Amplitude}} $$\n",
    "\n",
    "This can be plotted as a function of the base displacement frequency. A car suspension designer may want this ratio to be an optimal value for rider comfort. Maybe they'd like to make the ratio 1 or maybe even less than one if possible.\n",
    "\n",
    "**Exercise**\n",
    "\n",
    "Use the curve fitting technique from the previous notebook to plot $X/Y$ for a range of frequencies. Your code  should look something like:\n",
    "\n",
    "```python\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "def cosine_func(times, amp, freq, phase_angle):\n",
    "    return amp * np.cos(freq * times - phase_angle)\n",
    "\n",
    "frequencies = np.linspace(1.0, 20.0, num=100)\n",
    "   \n",
    "amplitudes = []\n",
    "    \n",
    "for omega in frequencies:\n",
    "    # your code here\n",
    "\n",
    "amplitudes = np.array(amplitudes)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, sharex=True)\n",
    "ax.set_xlabel('$\\omega_b$ [rad/s]')\n",
    "ax.set_ylabel('Displacement Transmissibility') \n",
    "\n",
    "ax.axvline(, color='black') # natural frequency\n",
    "ax.plot()#?\n",
    "ax.grid();\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "def cosine_func(times, amp, freq, phase_angle):\n",
    "    return amp * np.cos(freq * times - phase_angle)\n",
    "frequencies = np.linspace(1.0, 20.0, num=100)\n",
    "   \n",
    "amplitudes = []\n",
    "    \n",
    "for omega in frequencies:\n",
    "    traj = sys.sinusoidal_base_displacing_response(Y, omega, 20.0)\n",
    "    popt, pcov = curve_fit(cosine_func,\n",
    "                           traj[10:].index, traj[10:].car_vertical_position,\n",
    "                           p0=(Y, omega, 0.05))\n",
    "    amplitudes.append(abs(popt[0]))\n",
    "\n",
    "amplitudes = np.array(amplitudes)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, sharex=True)\n",
    "ax.set_xlabel('$\\omega_b$ [rad/s]')\n",
    "ax.set_ylabel('Displacement Transmissibility') \n",
    "\n",
    "ax.axvline(np.sqrt(sys.constants['suspension_stiffness'] / sys.constants['sprung_mass']), color='black')\n",
    "ax.plot(frequencies, amplitudes / Y)\n",
    "ax.grid();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# write you answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "The second thing to investigate is the *force transmissibility*. This is the ratio of the force applied by the suspension to the sprung car. Riders will feel this force when the car travels over bumps. Reducing this is also preferrable. The force applied to the car can be compared to the \n",
    "\n",
    "**Excersice**\n",
    "\n",
    "Create a measurement to calculate the force applied to the car by the suspension. Simulate the system with $Y=0.01$ m, $v = 10$ m/s, and the distance between bump peaks as $6$ m. Plot the trajectories.\n",
    "\n",
    "```python\n",
    "def force_on_car(suspension_damping, suspension_stiffness,\n",
    "                 car_vertical_position, car_vertical_velocity, travel_speed, time):\n",
    "    # write this code\n",
    "\n",
    "sys.add_measurement('force_on_car', force_on_car)\n",
    "\n",
    "# write code for Y and omega_b, etc\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "Y = 0.01  # m\n",
    "bump_distance = 6  # m\n",
    "\n",
    "def force_on_car(suspension_damping, suspension_stiffness,\n",
    "                 car_vertical_position, car_vertical_velocity,\n",
    "                 travel_speed, time):\n",
    "    \n",
    "    wb = travel_speed / bump_distance * 2 * np.pi\n",
    "    \n",
    "    y = Y * np.sin(wb * time)\n",
    "    yd = Y * wb * np.cos(wb * time)\n",
    "    \n",
    "    return (suspension_damping * (car_vertical_velocity - yd) +\n",
    "            suspension_stiffness * (car_vertical_position - y))\n",
    "\n",
    "sys.add_measurement('force_on_car', force_on_car)\n",
    "\n",
    "v = 10.0\n",
    "sys.constants['travel_speed'] = v\n",
    "wb = v / bump_distance * 2 * np.pi  # rad /s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = sys.sinusoidal_base_displacing_response(Y, wb, 10.0)\n",
    "traj[['car_vertical_position', 'car_vertical_velocity', 'force_on_car']].plot(subplots=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# write your answer here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.animate_configuration(fps=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Force transmissibility will be visited more in your next homework."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Arbitrary Periodic Forcing (Fourier Series)\n",
    "\n",
    "Fourier discovered that any periodic function with a period $T$ can be described by an infinite series of sums of sines and cosines. See the wikipedia article for more info (https://en.wikipedia.org/wiki/Fourier_series). The key equation is this:\n",
    "\n",
    "$$ F(t) = \\frac{a_0}{2} + \\sum_{n=1}^\\infty (a_n \\cos n\\omega_T t + b_n \\sin n \\omega_T t)$$\n",
    "\n",
    "The terms $a_0, a_n, b_n$ are called the Fourier coefficients and are defined as such:\n",
    "\n",
    "$$ a_0 = \\frac{2}{T} \\int_0^T F(t) dt$$\n",
    "\n",
    "$$ a_n = \\frac{2}{T} \\int_0^T F(t) \\cos n \\omega_T t dt \\quad \\textrm{for} \\quad n = 1, 2, \\ldots $$\n",
    "\n",
    "$$ b_n = \\frac{2}{T} \\int_0^T F(t) \\sin n \\omega_T t dt \\quad \\textrm{for} \\quad n = 1, 2, \\ldots $$\n",
    "\n",
    "\n",
    "## Introduction to SymPy\n",
    "\n",
    "SymPy is a Python package for symbolic computing. It can do many symbolic operations, for instance, integration, differentiation, linear algebra, etc. See http://sympy.org for more details of the features and the documentation. Today we will cover how to do integrals using SymPy and use it to find the Fourier series that represents a sawtooth function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import sympy as sm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `init_printing()` enables LaTeX based rendering in the Jupyter notebook of all SymPy objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sm.init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Symbols can be created by using the `symbols()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x, y, z = sm.symbols('x, y, z')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x, y, z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `integrate()` function allows you to do symbolic indefinite or definite integrals. Note that the constants of integration are not included in indefinite integrals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sm.integrate(x * y, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Integral` class creates and unevaluated integral, where as the `integrate()` function automatically evaluates the integral."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "expr = sm.Integral(x * y, x)\n",
    "expr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To evaluate the unevaluated form you call the `.doit()` method. Note that all unevaluated SymPy objects have this method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "expr.doit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This shows how to create an unevaluated definite integral, store it in a variable, and then evaluate it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "expr = sm.Integral(x * y, (x, 0, 5))\n",
    "expr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "expr.doit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fourier Coefficients for the Sawtooth function\n",
    "\n",
    "Now let's compute the Fourier coefficients for a saw tooth function. The function that describes the saw tooth is:\n",
    "\n",
    "$$\n",
    "F(t) = \n",
    "\\begin{cases} \n",
    "      A \\left( \\frac{4t}{T} - 1 \\right) & 0 \\leq t \\leq T/2 \\\\\n",
    "      A \\left( 3 - \\frac{4t}{t} \\right) & T/2 \\leq t \\leq T \n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "where:\n",
    "\n",
    "- $A$ is the amplitude of the saw tooth\n",
    "- $T$ is the period of the saw tooth\n",
    "- $\\omega_T$ is the frequency of the saw tooth, i.e. $\\omega_T = \\frac{2\\pi}{T}$\n",
    "- $t$ is time\n",
    "\n",
    "This is a piecewise function with two parts from $t=0$ to $t=T$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A, T, wT, t = sm.symbols('A, T, omega_T, t', real=True, positive=True)\n",
    "A, T, wT, t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first Fourier coefficient $a_0$ describes the average value of the periodic function. and is:\n",
    "\n",
    "$$a_0 = \\frac{2}{T} \\int_0^T F(t) dt$$\n",
    "\n",
    "This integral will have to be done in two parts:\n",
    "\n",
    "$$a_0 = a_{01} + a_{02} = \\frac{2}{T} \\int_0^{T/2} F(t) dt + \\frac{2}{T} \\int_{T/2}^T F(t) dt$$\n",
    "\n",
    "These two integrals are evaluated below. Note that $a_0$ evaluates to zero. This is because the average of our function is 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ao_1 = 2 / T * sm.Integral(A * (4 * t / T - 1), (t, 0, T / 2))\n",
    "ao_1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ao_1.doit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ao_2 = 2 / T * sm.Integral(A * (3 - 4 * t / T), (t, T / 2, T))\n",
    "ao_2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ao_2.doit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But SymPy can also handle piecewise directly. The following shows how to define a piecewise function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "F_1 = A * (4 * t / T - 1)\n",
    "F_2 = A * (3 - 4 * t / T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "F = sm.Piecewise((F_1, t<=T/2),\n",
    "                 (F_2, T/2<t))\n",
    "F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "F_of_t_only = F.xreplace({A: 0.01, T: 2 * sm.pi / wb})\n",
    "F_of_t_only"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sm.plot(F_of_t_only, (t, 0, 2 * np.pi / wb))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The integral can be taken of the entire piecewise function in one call."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sm.integrate(F, (t, 0, T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the Fourier coefficients $a_n$ and $b_n$ can be computed.\n",
    "\n",
    "$$\n",
    "a_n = \\frac{2}{T}\\int_0^T F(t) \\cos n\\omega_Tt dt \\\\\n",
    "b_n = \\frac{2}{T}\\int_0^T F(t) \\sin n\\omega_Tt dt\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = sm.symbols('n', real=True, positive=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For $a_n$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "an = 2 / T * sm.Integral(F * sm.cos(n * wT * t), (t, 0, T))\n",
    "an"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "an.doit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This can be simplified:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "an = an.doit().simplify()\n",
    "an"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now substitute the $2\\pi/T$ for $\\omega_T$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "an = an.subs({wT: 2 * sm.pi / T})\n",
    "an"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how this function varies with increasing $n$. We will use a loop but the SymPy expressions will not automatically display because they are inside a loop. So we need to use SymPy's `latex()` function and the IPython display tools. SymPy's `latex()` function  transforms the SymPy expression into a string of matching LaTeX commands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sm.latex(an, mode='inline')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `display()` and `LaTeX()` functions then turn the LaTeX string in to a displayed version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from IPython.display import display, Latex "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can see how $a_n$ varies with $n=1,2,\\ldots$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for n_i in range(1, 6):\n",
    "    ans = an.subs({n: n_i})\n",
    "    display(Latex('$a_{} = $'.format(n_i) + sm.latex(ans, mode='inline')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For even $n$ values the coefficient is zero and for even values it varies with the inverse of $n^2$. More precisely:\n",
    "\n",
    "$$\n",
    "a_n =\n",
    "\\begin{cases}\n",
    "0 & \\textrm{if }n\\textrm{ is even} \\\\\n",
    "-\\frac{8A}{n^2\\pi^2} & \\textrm{if }n\\textrm{ is odd}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "SymPy can actually reduce this further if your set the assumption that $n$ is an integer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = sm.symbols('n', real=True, positive=True, integer=True)\n",
    "an = 2 / T * sm.Integral(F * sm.cos(n * wT * t), (t, 0, T))\n",
    "an = an.doit().simplify()\n",
    "an.subs({wT: 2 * sm.pi / T})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The odd and even versions can be computed by setting the respective assumptions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = sm.symbols('n', real=True, positive=True, integer=True, odd=True)\n",
    "an = 2 / T * sm.Integral(F * sm.cos(n * wT * t), (t, 0, T))\n",
    "an = an.doit().simplify()\n",
    "an.subs({wT: 2 * sm.pi / T})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that $b_n$ is always zero:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bn = 2 / T * sm.Integral(F * sm.sin(n * wT * t), (t, 0, T))\n",
    "bn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bn.doit().simplify().subs({wT: 2 * sm.pi / T})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical evalution of the Fourier Series\n",
    "\n",
    "Now the Fourier coefficients can be used to plot the approximation of the saw tooth forcing function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function plots the actual sawtooth function. It does it all in one line by cleverly using the absolute value and the modulo functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sawtooth(A, T, t):\n",
    "    return (4 * A / T) * (T / 2 - np.abs(t % T - T / 2) ) - A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = 1\n",
    "T = 2\n",
    "t = np.linspace(0, 5, num=500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "\n",
    "plt.plot(t, sawtooth(A, T, t));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "# Exercise\n",
    "\n",
    "Write a function that computes the Fourier approximation of the sawtooth function for a given value of $n$, i.e. using a finite number of terms. Then plot it for $n=2, 4, 6, 8, 10$ on top of the actual sawtooth function. How many terms of the infinite series are needed to get a good sawtooth?\n",
    "\n",
    "```python\n",
    "def sawtooth_approximation(n, A, T, t):\n",
    "    # code here\n",
    "    return f\n",
    "\n",
    "# plot sawtooth\n",
    "f = sawtooth(A, T, t)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(t, f, color='k', label='true sawtooth')\n",
    "\n",
    "for n in np.arange(2, 12, 2):\n",
    "    f_approx = sawtooth_approximation(n, A, T, t)\n",
    "    plt.plot(t, f_approx, label='n = {}'.format(n))\n",
    "\n",
    "plt.legend()\n",
    "# zoom in a bit on the interesting bit\n",
    "plt.xlim(0, T)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "def sawtooth_approximation(n, A, T, t):\n",
    "    # odd values of indexing variable up to n\n",
    "    n = np.arange(1, n+1)[:, np.newaxis]\n",
    "    # cos coefficients\n",
    "    an = A *(8 * (-1)**n - 8) / 2 / np.pi**2 / n**2\n",
    "    # sawtooth frequency\n",
    "    wT = 2 * np.pi / T\n",
    "    # sum of n cos functions\n",
    "    f = np.sum(an * np.cos(n * wT * t), axis=0)\n",
    "    return f\n",
    "\n",
    "# plot sawtooth\n",
    "f = sawtooth(A, T, t)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(t, f, color='k', label='true sawtooth')\n",
    "\n",
    "for n in np.arange(2, 12, 2):\n",
    "    f_approx = sawtooth_approximation(n, A, T, t)\n",
    "    plt.plot(t, f_approx, label='n = {}'.format(n))\n",
    "\n",
    "plt.legend()\n",
    "# zoom in a bit on the interesting bit\n",
    "plt.xlim(0, T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# write answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is a interactive plot that shows the same thing as above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = 1\n",
    "T = 2\n",
    "t = np.linspace(0, 5, num=500)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "\n",
    "f = sawtooth(A, T, t)\n",
    "\n",
    "saw_tooth_lines = ax.plot(t, f, color='k')\n",
    "\n",
    "n = 2\n",
    "\n",
    "f_approx = sawtooth_approximation(n, A, T, t)\n",
    "\n",
    "approx_lines = ax.plot(t, f_approx)\n",
    "\n",
    "leg = ax.legend(['true', 'approx, n = {}'.format(n)])\n",
    "\n",
    "# zoom in a bit on the interesting bit\n",
    "plt.xlim(0, 2 * T)\n",
    "\n",
    "def update(n=0):\n",
    "    f_approx = sawtooth_approximation(n, A, T, t)\n",
    "    approx_lines[0].set_ydata(f_approx)\n",
    "    leg.get_texts()[1].set_text('approx, n = {}'.format(n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from ipywidgets import interact\n",
    "\n",
    "interact(update, n=(0, 20, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "# Apply the sawtooth to the quarter car\n",
    "\n",
    "Now that you know the Fourier series coefficients. Calculate them for a suitable number of terms and simulate them with the `sys.periodic_base_displacing_response()` function.\n",
    "\n",
    "Your code should look something like:\n",
    "\n",
    "```python\n",
    "def fourier_coeffs(A, T, N):\n",
    "    # write your code here\n",
    "\n",
    "a0, an, bn = fourier_coeffs(?)\n",
    "\n",
    "traj = sys.periodic_base_displacing_response(?)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "def fourier_coeffs(A, T, N):\n",
    "    n = np.arange(1, N+1)\n",
    "    an = A *(8 * (-1)**n - 8) / 2 / np.pi**2 / n**2\n",
    "    return 0, an, np.zeros_like(an)\n",
    "\n",
    "a0, an, bn = fourier_coeffs(0.01, 2 * np.pi / wb, 100)\n",
    "\n",
    "traj = sys.periodic_base_displacing_response(a0, an, bn, wb, 20.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "traj.plot(subplots=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sys.animate_configuration(fps=30)"
   ]
  }
 ],
 "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
}