Semi_active Damping.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import the SymPy Module\n",
    "import sympy\n",
    "\n",
    "# Import the necessary sub-modules and methods for dynamics\n",
    "from sympy.physics.mechanics import dynamicsymbols\n",
    "from sympy.physics.mechanics import inertia, LagrangesMethod, Lagrangian\n",
    "from sympy.physics.mechanics import Particle, Point, ReferenceFrame, RigidBody\n",
    "\n",
    "# initiate better printing of SymPy results\n",
    "sympy.init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "# Define the genearlized coordinate\n",
    "x, theta1, f = dynamicsymbols('x theta_1 f')\n",
    "\n",
    "# Also define the first derivative\n",
    "x_dot, theta1_dot = dynamicsymbols('x theta_1', 1)\n",
    "\n",
    "# Define the symbols for the other paramters\n",
    "m, mt, g, k, L, Izz, t = sympy.symbols('m mt g k L I_zz t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N = ReferenceFrame('N')\n",
    "\n",
    "# Define a body-fixed frame along the spring, with y aligned from m to the pin\n",
    "A = N.orientnew('A', 'Axis', [theta1, N.z])\n",
    "# B = N.orientnew('B', 'Axis', [theta2, N.z])\n",
    "\n",
    "# Define an origin and set its velocity to 0\n",
    "O = Point('0')\n",
    "O.set_vel(N, 0 * N.x)\n",
    "\n",
    "# Create the trolley - treat as a particle\n",
    "T = O.locatenew('T', x * N.x)\n",
    "T.set_vel(N, x_dot * N.x)\n",
    "\n",
    "# Define the point at the end of the spring and set its velocity\n",
    "P = T.locatenew('P', -L * A.y)\n",
    "P.set_vel(N, x_dot * N.x + L * theta1_dot * A.x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAscAAAAeCAYAAADAWLPgAAAABHNCSVQICAgIfAhkiAAADyhJREFU\neJztnXmwHUUVh7+QnQBhC1sMZCFIQAMBEjAk5IWwSBCssAZBLCFgCQKJoBQoiKAgCIqAypLIk0VA\nUiIEBEEgmqDsWBolaLlgMBAgBgFZQiD+8Zup23fu3Lmz9J2ZO6+/qlfv3XlzZzm/0326e073gMPh\ncDgcDofD4agcs4HXgVFFX4jD4XBUjCrUr1W4h7LgbFltnL6W2B2YA5wHPADsVcA1jAL+WcB505DV\nXofYviCHw+GIoJPq12ZU4R7KQhlt6eKqPaqmbyHargdcZHw+AngLGJrzdZwCXJvzOdOQ1V7bAz+3\nfVEOh8MRQafUr1FU4R7KQtls6eKqXaqkrzVt1wHuAp4F1gJrgEdRS/1jIfuPBT6gNvy+gfe9I2xc\nTAR9gAuAE4GTgUeAw2J+dzhwDXAxMB9YP8X5twIGtNhnM+CzgW1Z7XU+cGRg2wTgZ8C9wB+BeeTf\nOUnqNz2BKcDEkO3DKc7/+gGXAn1TnNMGPclPOln/IuvXtPcWJMs9JMHW9ZaZstuyqnE1L/LQN0s5\nyaKvdW27vJN/q8V+vdBwdy/v847e98bFPVFK5qK8GICNUZDdMMb3hgNPA1t6n08Hrkh47sHAqSHb\nRwHLgW2MbTOAScbnrPZ6GhhofN4FuJ/ava8H/AZ4Gd1rEmw0mLqI5zdVZzRqHAQZTrH+B/K1SxKe\n08T5SWvy1j9Me0ivf1H1a1bfNkl7D0nIqkUciurImpTdllWPq1lpZkefduubtVxn0de6tud4J58W\nZ2eDG4HLEn4nKWOBN6nd8N7A72J8rx/wFHCcse1w4F8Jz/91ZNAgs4F3qO8d9QKuizhWEntN8PY3\nuQfYNrBtHNLu1pjH9bmO7D3jtH5TNW5HFYJJmfzvSlTZpMH5SWvy1j9Me0inf5H1q626Ne09JMWG\nFq2wUd6y0Im2rFpczUozO0I++tpsM0F8fVNru07EQacA75HMSMcBLwJnJPhOGvYBFgNve5+noZ7A\nRi2+dxoa0bjJ2DYYGEb83t0QoD9ypiCTgceQ2D5rgVcJn/WZ1F5HAz8JbJsCPIQeR/g8A7yG7JSE\nvmTv5abxm6oxEXgF+G9ge5n870rgzJjnDOL8JJoi9A/THtLpX1T9arNuTXsPSbClRStslLcsdJot\nqxhXs9LMjtB+fW2Wa0imb2ptmzWO+6L8v6dQ0nMcpnu/v4wMMTzm99KwChkH1Bs5FHgQGaIZA1Aw\nmAesNraP8X5HdRRM9kf5OGFMQsPzQR5Fj5FNktqrNzAV5Waa/APYHBgU2P4u9Y8S8iCN31SRk2js\nrZbN//4CjMBugItL1f2kCP2baQ/J9S+qfrXl25DuHpJiS4uy00m2rGJctUGUHdutr81ynUTfTNr2\naXLQ3YB1aW7MIHuhEYN7gC2APYAV1C8F8knUKt8J+AzKazkc9RImokki9wJfBDZBvY1+KDn7vcD5\nbgH2BI7ybuZmZMinIq5xpnfc4CORPYE3kGHiMJH6XMIjvWvcCPVEpgL3AQuA73v7LKc+hyaOvYLs\nDSxCuUAme6AJLy8Z27ZC4i+MdUf2SOo3II3PRjmarwBLgX2BO9HImy2f2Rk4F/UO30HBei7wpLHP\n7miW7gpUNjYELgSeM/aZgRp2q1Dh3BQ9BjrZ2GcS9Y+WoVz+5/M4sB9wW8xz2yKun9isM6qofxzt\nIbn+RdWvNn07zT0A7Aocg3xsODAL+BzyhaHA14C/hVxvVi1sU3S8hfxtWaW4akO/uHZst762ynVS\nfdui7Znohg9stSMwEi0WvTbws4GxTz/gcu/vJ9AQ/unUkqvPRMnQ3/GOB2r1v4EcwwZ3oscG9xk/\nDwLvIwPG5Q4aexwgp2jWqxxNLYcmjr3CuB45cBwuRPcVNlM+im6yjfgn8RuQrRZSH5zPQs58KvZ8\nZiKyeZex7QbgeePzQagRNMTYNgYVvHHG54cCx56J7OYzmvoGl09Z/M/kJOqXx4lLN+33E5t1RtX1\nj9Levyab+odh4/5s+nYaRqN0E9/HutEI+0RU936AfNDHphZRdJOsvJUh3uZtyyrFVdv6tbJjGpLo\na6Ncp9E3k7bNRo6noJtbHOOgf29xgf7xFiNDjgQepj6Z+j3U+7nFOx7ehb6PWvNZ6e1dwx3Ap4zt\n01Hv4uEExxpE+CjIVOTIb4f8byByXIhnryADUID+bYx9R6LRr4ti7m+TJH4D6k3uABxsbHsG6fUu\ndnymFyoki6jvFa5EuU6gR0nXA99Ao9c+z6Kgfz0aeRyLerrrU9NzATDe+M5QVHGZlMn/TP5DMY94\n4/iJrTqjJ+gfpT20X39b92fTt9MwB/gSCrr+9axC9egw1BjpNva3qYVNio63kL8tqxRXbevXyo5p\nSKKvjXKdVN/M2oY1jnuj1vYfaJxQ4vNRtDZcXJagx5lj0eOBywP/H48e8T0RuODBwJ8SnKcZQ71j\nBSf/HOD9vj3BsVZ6x1oZ2N6F8u3C2JhaTk8aDkSPUNa22K8/Sj6/DvhqxH4/Ro9rgmyNZneuDvnf\n8UQ/ZknqN8OAz6PHKK8b+0xCtroLOz6zO7AdjY+O5xh/H4gqm6Uh1/wcGsXeFTWwNgNe8K7vcfQY\n2RwBGULj/ZfV/1bRuKKCSZF+YqvO6An6d9Fce0ivf1xs3V8RdavJxcD/jM8TqQX4ZTROALKpBdgr\nb0XHWyjelnEoa1y1rV8X0XZMQxJ9q9BmAhQE1gLfa/L/gcBV8a+xjtnIoP0C25ejYW2Tk1CPImz5\nj6SMp3HZqD4oXyWY97g+Wrh+6ybHugAFbxN/zb39mnxnDtlecDAf+EiM/W5EC16npZv0j8uT+s1s\nb//9A/ston52aVafmemd58SIaz/L2ydsFvIJ3v/8RcRHoBccLPO2v4NytXyORD16kyT+B9E+aNP/\npqORz6R0k6+f9CT9k9Y/rbQH+/oHsXV/RdStzRhDc5/wsalFFN2kK29FxdsgZbKlSdnjqg394tgx\nK6307cg2U9gMYv991b9ucrDjUT5ZGqaiWYtmD2o7lGS9MLDvDO88b6JeUTDXJMmPn5BtJmAfgEaB\nzjG2zUKjQIfSfHb1QzSKNtU7hz8kvyEaGfUZQ30vLwmDUSBZ0mK/81DOz7nGtmNTnjMNSf1me+/3\nY8a2ASjQLjS2ZfEZUEUC0asy+D3UzUL+t4n3ewXqyb+LcqWGIV1uQ42l/t5+r9K4eHpc/4PWPmjT\n/zaiPo0gD5L6SU/SP03900p7iKd/2rrV5v1l9e0sMSLINORz5mPZkYF9bGrRDoqKt0HKaMtOiKtZ\n6z7/GGF2tKUttNa3I9tMYRXUFO932CSKzdFISprGcW8UGIP5Z10oh8Zc6mNjZLybvc9+MnraH38R\nenPW4uloKN0M0nORwaJYhIKkyWTg99TW8TvNONcIJFJwxmRcDkWvOYziGKTlBYHtecyK9knqN6+h\nnKLXAsfoT81HsvoMqPC9gHIfg8xAs14XoCXFPhyyz67e9xejN+vMNP63DK1K8AG1BtGL1BpUPnH9\nD1r7oE3/25Ra4zEvkvhJT9M/Tf0TpT3E1z9t3Wrz/rL6dpYYMRC9NdAf4doXpf74Sw2uQ2MqgE0t\nbFNkvO0EW5Y9rtqo+6C5HbOU96T6dmSbKdg47oMueimNI0qTUYAIWxojDuNQAFkY2N6Feghm/spw\n5Bz3o2D6eIrzmaxEQdIfrTwONcLCXmfYijXUAqVPb2rLiYxHjzdeRM50BHB1ivP4zKTxMa3JXsB3\nUU/tJuPnVvQIMw/S+M18tJyXP1q3DfBD4N/AX71tNnxmDdJ7Cloex2cIKtQvIf84EY1abmnsMwJV\nOsca134G9SsafAj1Pld4n5eiRp6Zz19W/9uN8JUV2kVSP3H6NxLUv5n2kJ/+tu4v77rVZDrSdkd0\nHyOpn0T0FRrXrraphW2KjLedYMuyx1Vb+kXZMS1J9e3oNtP2wK9Qo2QtmlDySzSC8zBaNNkfVp+Q\n8oIPRr2L4CTAB9DELJPeaBLH1ShJuhfhzEaTucw3qWxHeK7KDsDd3jGvIPw1ij5raZ0jdLxxjJ1Q\nT+5SNIPT73RsQrb3k29B88fPPitp/ggk2CtqRTfJcqOy+s0JSJNL0IjSk9S/Xcumz0wAfoHymS9D\nOVvByUiTvf//AC1TM4/6XKljUMG/DPgmWt3gKhof9/0UjTiaJPE/aO2DNvxvCekq+u4W1xYkrZ/0\nVP2T1D/NtIfs+gfr12Z1K9i9P5t1a9wYsSnwI+Db3s8gNMHqGuQL+8a43ixaRNFN8pzVIuNtmW0J\n5Y+rYE+/KDsGiVve0+pbxTZTaRlF4wLQ81CjKwtxglMezEaJ9nnRTXH3vS7KW5pV0Plt8nHUeMpC\nu31wPBqpT0M35SgfZSWr/nnUP3H0D9avNupWyLd+bVeMyJNuylHeqmBL6FlxNQntKu95kre2peUU\n4NrAtufRI/oslKVx/AjqpeXFLuTzDvgh1K9tDHAYWqcxOArXqdxNtkXX2+2DN5C+nOTlJ51MFv3z\nqH/i6B+sX23UrZBv/dquGJEnZSlvVbAlVDeuZqVd5T1PrGnb7CUgZaUPeiXhMuRsM1HuCCgP5gA0\nPH8UynGJ+zKKMrItmqz2ao7nfDqn81wJHIIen7yF8jwvQY+EluV0De3mfDRT/+yiLySE/VFu7PMp\nv5+Xn3Qynah/s/q1k+rWKsaIospbFW1Z5bialCqUd5MitC0Nc9GwOWiG5hrq81NmkS2R+2j0qHEt\nSs7+QoZjZeVclONYRY5C+Z8XoZyj+cAnCr2i9jCT5G8ga7cP9kXr+jbLK3TYI6n+edQ/UfpH1a9Z\n61bI5/7aHSN6ElW0ZZXjalLaXd7zpsdqOxYt++E/qtybxjcy3UBtsf5OZwHtWZDd4XA4grSqXzuh\nbu1pMaKdVNWWLq6KKpT3IFa1jZq5WDb2QUP7/nu4p6GlS8zF/Seh2e/r0LjWaKdxELU1AB0Oh6Od\ntKpfO6Fu7Wkxop1U1ZYurooqlPcgVrXtpMbxKmrr862HFnp+ED2qAwm4GngZ+DTZJkM5HA5HTyKq\nfu2UutXFCHs4W1abKpT3ttK76AtIwFKUJD4I2BktOTIC+DNaP3U1MBGtXboceKaQq3Q4HI7OI6p+\nXUJn1K0uRtjD2bLaVKG8OxwOh8PhcDgcDofD4XA4HA6Hw+FwOBwOh8PhcDgcDofD4XA4HA5HYfwf\nDn+Ar+S0PgwAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\frac{L^{2} m}{2} \\frac{d}{d t} \\theta_{1}{\\left (t \\right )}^{2} + L g m \\cos{\\left (\\theta_{1}{\\left (t \\right )} \\right )} + L m \\cos{\\left (\\theta_{1}{\\left (t \\right )} \\right )} \\frac{d}{d t} \\theta_{1}{\\left (t \\right )} \\frac{d}{d t} x{\\left (t \\right )} + \\frac{m}{2} \\frac{d}{d t} x{\\left (t \\right )}^{2} + \\frac{mt}{2} \\frac{d}{d t} x{\\left (t \\right )}^{2}$$"
      ],
      "text/plain": [
       "                2                                                             \n",
       " 2   ⎛d        ⎞                                                             ⎛\n",
       "L ⋅m⋅⎜──(θ₁(t))⎟                                                           m⋅⎜\n",
       "     ⎝dt       ⎠                                      d         d            ⎝\n",
       "───────────────── + L⋅g⋅m⋅cos(θ₁(t)) + L⋅m⋅cos(θ₁(t))⋅──(θ₁(t))⋅──(x(t)) + ───\n",
       "        2                                             dt        dt            \n",
       "\n",
       "         2                2\n",
       "d       ⎞       ⎛d       ⎞ \n",
       "──(x(t))⎟    mt⋅⎜──(x(t))⎟ \n",
       "dt      ⎠       ⎝dt      ⎠ \n",
       "────────── + ──────────────\n",
       "   2               2       "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Define the COM of the bar location and velocity\n",
    "# G = P.locatenew('G', -l/2 * B.y)\n",
    "# G.v2pt_theory(P, N, A)\n",
    "\n",
    "# Create the bar rigid bogy\n",
    "# I_bar = inertia(B, 0, 0, Izz)\n",
    "# bar = RigidBody('bar', G, B, m, (I_bar, G))\n",
    "\n",
    "# Treat the trolley as a particle\n",
    "mtr = Particle('mtr', T, mt)\n",
    "end = Particle('end', P, m)\n",
    "\n",
    "# Now, set the potential energy of the bar\n",
    "h = P.pos_from(T) & N.y # Define the height of the COM by the vertical distance from the datum\n",
    "end.potential_energy = (m * g * h)\n",
    "\n",
    "# # Set up the force list - each item follows the form:\n",
    "# #    (the location where the force is applied, its magnitude and direction)\n",
    "# # Here, there are no non-conservataive external forces\n",
    "forces = [(T, f * N.x)]\n",
    "\n",
    "# # Form the Lagrangian\n",
    "L =  Lagrangian(N, mtr, end)\n",
    "\n",
    "# # Print the Lagrangian as a check\n",
    "L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "LM = LagrangesMethod(L, [x, theta1], forcelist = forces, frame = N)\n",
    "EqMotion = LM.form_lagranges_equations()\n",
    "lrhs = LM.rhs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# import NumPy with namespace np\n",
    "import numpy as np\n",
    "\n",
    "# import the scipy ODE solver\n",
    "from scipy.integrate import ode\n",
    "\n",
    "# import the plotting functions from matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# set up the notebook to display the plots inline\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define the forcing function\n",
    "def force(t):\n",
    "    # set the maximum force output\n",
    "    fmax = 0.0\n",
    "    \n",
    "    # return a bang-bang command\n",
    "    return fmax * (t > 0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define the states and state vector\n",
    "w1, w2, w3, w4 = sympy.symbols('w1 w2 w3 w4', cls=sympy.Function)\n",
    "w = [w1(t), w2(t), w3(t), w4(t)]\n",
    "\n",
    "# Set up the state definitions and parameter substitution\n",
    "sub_params = {x : w1(t),\n",
    "              theta1: w2(t),\n",
    "              x_dot : w3(t),\n",
    "              theta1_dot : w4(t), \n",
    "              m : 10.0,\n",
    "              mt : 100.0,\n",
    "              g : 9.81,\n",
    "              L : 2.0,\n",
    "              f: force(t)}\n",
    "\n",
    "# set this parameter to enable array output from sympy.lambdify\n",
    "mat2array = [{'ImmutableMatrix': np.array}, 'numpy']\n",
    "\n",
    "# Create a function from the equations of motion\n",
    "# Here, we substitude the states and parameters as appropriate prior to the lamdification\n",
    "eq_of_motion = sympy.lambdify((t, w), \n",
    "                              lrhs.subs(sub_params), \n",
    "                              modules = mat2array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Set up the initial conditions for the solver\n",
    "x_init = 0.0                   # Initial Trolley position (m)\n",
    "x_dot_init = 0.0                # Initial trolley velocity (m/s)\n",
    "\n",
    "theta1_init = 0.7               # Initial angle (rad)\n",
    "theta1_dot_init = 0.0           # Initial angular velocity (rad/s)\n",
    "\n",
    "\n",
    "# Pack the initial conditions into an array\n",
    "x0 = [x_init, theta1_init, x_dot_init, theta1_dot_init]\n",
    "\n",
    "# Create the time samples for the output of the ODE solver\n",
    "sim_time = np.linspace(0.0, 30.0, 900) # 0-10s with 1001 points in between"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "can't convert expression to float",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-21-666d4dbccbe2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     16\u001b[0m \u001b[0;32mwhile\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msuccessful\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mt\u001b[0m \u001b[0;34m<=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0msim_time\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m1e-9\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     17\u001b[0m     \u001b[0mresponse\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 18\u001b[0;31m     \u001b[0mr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mintegrate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mt\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mdt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     19\u001b[0m     \u001b[0mindex\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m//anaconda/lib/python3.5/site-packages/scipy/integrate/_ode.py\u001b[0m in \u001b[0;36mintegrate\u001b[0;34m(self, t, step, relax)\u001b[0m\n\u001b[1;32m    406\u001b[0m             self._y, self.t = mth(self.f, self.jac or (lambda: None),\n\u001b[1;32m    407\u001b[0m                                 \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_y\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 408\u001b[0;31m                                 self.f_params, self.jac_params)\n\u001b[0m\u001b[1;32m    409\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mSystemError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    410\u001b[0m             \u001b[0;31m# f2py issue with tuple returns, see ticket 1187.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m//anaconda/lib/python3.5/site-packages/scipy/integrate/_ode.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, f, jac, y0, t0, t1, f_params, jac_params)\u001b[0m\n\u001b[1;32m    863\u001b[0m         args = ((f, jac, y0, t0, t1) + tuple(self.call_args) +\n\u001b[1;32m    864\u001b[0m                 (f_params, jac_params))\n\u001b[0;32m--> 865\u001b[0;31m         \u001b[0my1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mistate\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrunner\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    866\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mistate\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    867\u001b[0m             warnings.warn(self.__class__.__name__ + ': ' +\n",
      "\u001b[0;32m//anaconda/lib/python3.5/site-packages/sympy-0.7.7.dev0-py3.5.egg/sympy/core/expr.py\u001b[0m in \u001b[0;36m__float__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    224\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_number\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mas_real_imag\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    225\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"can't convert complex to float\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 226\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"can't convert expression to float\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    227\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    228\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__complex__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: can't convert expression to float"
     ]
    }
   ],
   "source": [
    "# Set up the initial point for the ode solver\n",
    "r = ode(eq_of_motion).set_initial_value(x0, sim_time[0])\n",
    " \n",
    "# define the sample time\n",
    "dt = sim_time[1] - sim_time[0]   \n",
    "\n",
    "# pre-populate the response array with zeros\n",
    "response = np.zeros((len(sim_time), len(x0)))\n",
    "\n",
    "# Set the initial index to 0\n",
    "index = 0\n",
    "\n",
    "# Now, numerically integrate the ODE while:\n",
    "#   1. the last step was successful\n",
    "#   2. the current time is less than the desired simluation end time\n",
    "while r.successful() and r.t <= (sim_time[-1] + 1e-9):\n",
    "    response[index, :] = r.y\n",
    "    r.integrate(r.t + dt)\n",
    "    index += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAGVCAYAAACy6/bkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm8HGWV//HPyR5IAklIQiSEsO8RCGiAsAYQgRAQRPip\nqKj8HEeYAUeiMoji4A8BGXDc0DGIOCOCiAkSFiObEmSXfScJS0JI4JKNQJZ7fn9UddLpdN/b1V1P\nV1f39/169atvVVdXnxya5Nx6njqPuTsiIiIiko0eWQcgIiIi0s5UjImIiIhkSMWYiIiISIZUjImI\niIhkSMWYiIiISIZUjImIiIhkqFfWATSCmc0B+pTuBla5++jGRyQiIiISaYtiDOh09w9kHYSIiIhI\nKQ1TioiIiGRIxZiIiIhIhlSMiYiIiGSoXYoxM7MLzexJM5tvZk+b2Q/MbGjWgYmIiEh7a5dirBN4\nF9gXGAV8Gfg48KCZDc8yMBEREWlv5u5ZxxCcmQ1x97dL9h0DTAf+291PzyYyERERaXdtUYyVY2Y9\ngPeARWp7ISIiIllplz5jG3D3TjN7CxhW7nUz83POOYf+/fsDMG7cOCZMmMDgwYMB6OjoANB2gu2X\nX+7B3ntvwqBBnTz//GL69Gmu+Fphu7CvWeJppe1ly5ax5ZZbNk08rbT96quvMmDAgKaJp5W2S/9u\nyDqeFts2UtLyV8bM7CCgt7vPLNlfuDL2trtvXuZ93uq5abRLL4WvfQ1OPPF9rr++b9bhtKSOjo61\nf2FIupTbcJTbcJTboFIrxtphAv/BwBll9h9JdGXwloZG08ZuvDF6PukkFWKh6C/dcJTbcJTbcJTb\nfGiHYgzgGDP7spn1BjCzfYEfAfOB8zKNrE288Qbcdx/07QtHHpl1NCIiIs2jHYqxHwFnAycDs+N5\nYr8FbgX2dvfXsgyuXUyfDu5w2GGwenVH92+QmhTPD5F0KbfhKLfhKLf50PIT+N39LeCK+CEZ+eMf\no+fjj882DhERkWbT8hP4a6UJ/OlZsgSGDYPVq2H+fBiuNrsiknPXXnst5557LnPmzKGzszPrcCSA\nXr16MXHiRKZNm0bfvmXnOmsCv+THLbfAypWw//4qxEQk/6699lrOOussrrnmGlasWIG769GCjyVL\nltC3b18uuOCC4N8pFWMSXGGI8rjjomfNYQhHuQ1HuQ0nb7k999xzueGGG9hvv/3o06dP1uFIIP37\n9+eyyy5j6tSpwT9Lw5QVaJgyHe+/Hw1RLl0KL70E22yjvjchKbfhKLfh5C23PXv2ZMWKFSrE2sDq\n1avp27cva9asKfeymr6GpmIsHbfeCh/9KIwdC489lnU0IiL1MzP070P76OK/t+aMST4UGr0WhihF\nRERkfSrGJJjOTpg2Lfq5uKVF3uaH5IlyG45yG45yK1mYO3cuX/rSl/j617/OSSedxPLlyzOLRcWY\nBPP3v8OCBbDVVvDBD2YdjYiISGTOnDl87GMf44ILLuCiiy5i/PjxfPOb38wsHhVjEkxxo1crGlnP\n00TdvFFuw1Fuw1FupZFWrVrFiSeeyJlnnsnwuN/S6NGjmVYYysmAijEJwl3zxUREpPlcfvnlvPHG\nG3zyk59cu2/x4sW8+uqrle6aDE7FmATx9NPw4oswdGjU7LWY5oeEo9yGo9yGo9xKo7z//vtcfPHF\nfOELX6BXr3UrQj7zzDMAma2m0PJrU0o2CkOUxx4LvfQtExFpuEceeYTf/OY3mBlz587lF7/4BVde\neSXvvPMOr7/+OhdccAFbb7111mE21G9/+1vefvttPvGJT6y3/95772XgwIH07t07k7j0z6QE0dUQ\npeaHhKPchqPchqPcpu/FF1/k6quv5oorrgDgc5/7HOPHj+fqq6+ms7OTAw44gL322ouzzjor40gb\n649//CP9+vXjq1/9Ku6OmbFy5UoefPBB9i8dxmkgFWOSuldfhYcfho02gsMPzzoaEZH2c/nll3PJ\nJZes3V6+fDlDhgxh/PjxvPbaa3z1q1/ls5/9bHYBZqCzs5O7776bj33sY1xzzTVr999yyy3ccccd\nHHrooZnFpjljkrrCEOWRR0L//hu+rvkh4Si34Si34Si36ZsyZQr9i/4CnjVrFocddhgAo0aN4uKL\nL267K5Kvv/46ixcvZvz48evtnzFjBmbGiSeemFFkujImARS3tBARaTeW0iI59ay4tOWWW679+dln\nn2XevHkccsghKURVH/tO+eT4+eX/sKXHVzquGgsWLABgl112WbtvzZo1XH/99Rx44IHsuuuuNZ+7\nXirGJFVvvQV33x1N2j/66PLHtNtvY42k3Iaj3Iaj3IY1c+ZM+vbty3777bd23+zZs9tu8n6vXr0w\nMzbffPO1+2bMmMHChQv5/e9/n2FkKsYkZTffDGvWwGGHgf5+FZF2lPUa4u+99x7nn38+p556Krvu\nuiszZ85k7Nix9OvXL47PufTSS/nxj3/c8NiSXtmq50pYqdGjRwOs19Lisssu4/TTT2fChAmpfU4t\nNGdMUlUYouyq0avmh4Sj3Iaj3Iaj3KZrxowZXHrppTz11FM899xzvPzyy/Tt23ft6xdeeCGnnnpq\nhhFmY8iQIey33348++yzAEydOpWVK1euveM0S+ZZl/BNysxcuUnm3Xdhs81gxYrojspRo8of19HR\noWGJQJTbcJTbcPKWWzOjmf99eOutt5gyZQpDhw4F4Nvf/jZf/vKX6devH3369OHYY49l4sSJGUeZ\njWeeeYYpU6YwatQo+vTpw8UXX0yfPn26fE8X/71Tmh2oYqwiFWPJTZsWXRHbZx944IGsoxERCaPZ\nizFJVyOKMQ1TSmq0FqWIiEhyKsYkFatXw003RT9319JC80PCUW7DUW7DUW6l3akYk1T89a/w9tuw\nww6w005ZRyMiIpIfKsYkFcWNXrtreJinibp5o9yGo9yGo9xKu1MxJnVzr66lhYiIiGxIxZjU7dFH\n4ZVXYORI+NCHuj9e80PCUW7DUW7DyVtue/TowcqVK7MOQxpg9erV9OgRvlRSMSZ1K1wVmzwZGvCd\nFRHJ1JgxY3jooYeyDkMaYO7cuQwfPjz45+ifTqlb0pYWmh8SjnIbjnIbTt5ye+GFF3LCCScwa9Ys\nXSFrYStWrODss8/mtNNOC/5ZavpagZq+VufFF2H77WHQIFi4ELppZCwi0hKuvfZazj33XObMmUNn\nZ2fW4UgAvXr1YuLEiUybNm295aSKqAN/aCrGqnPppfC1r8Epp8D//m9178nb0id5otyGo9yGo9yG\no9wGpQ780hyKW1qIiIhIcroyVoGujHVvwYLoDso+faIhyoEDs45IRESkYXRlTLI3fXrUY+yww1SI\niYiI1ErFmNSs1kaveesplCfKbTjKbTjKbTjKbT6oGJOaLFkCM2dGSx9NmpR1NCIiIvmlYkxqcuut\nsHIl7L8/jBiR7L26sycc5TYc5TYc5TYc5TYfVIxJTZI2ehUREZHyVIxJYu+/DzffHP1cSzGmOQzh\nKLfhKLfhKLfhKLf5oGJMErvzTli6FHbfHbbdNutoRERE8k3FmCRWb6NXzWEIR7kNR7kNR7kNR7nN\nh7YrxszsA2a22MzWZB1LHnV2wrRp0c+aLyYiIlK/tivGgJ8CalFao/vvhzfegK22gj32qO0cmsMQ\njnIbjnIbjnIbjnKbD21VjJnZx4FdgQezjiWvihu9WmoLQYiIiLSvtlmb0sw2BZ4GTgXOBQ50955d\nHK+1KUu4w447wgsvRJP4Dz4464hEREQyo7Upa3ApcLu7z8w6kLx65pmoEBs6FCZMyDoaERGR1tAW\nxZiZHQwcA5yVcSi5Vmj0OmkS9OpV+3k0hyEc5TYc5TYc5TYc5TYfWr4YM7O+wJXAV91d38o61NvS\nQkRERDbU8nPGzOx7wF7ufmTRvjvRnLFEXn0VRo+GjTaCRYugf/+sIxIREclUanPG6hhsan5mNhb4\nMlBjEwYpKPQWO/JIFWIiIiJpauliDDgKcGCWrevDYMAQADObF++71N0vK33zlClT6B9XHuPGjWPC\nhAlruxkXxuHbZfv3v18F9F7b6LWe8xXPYWiWP1+rbBf2NUs8rbS9bNkyttxyy6aJp5W2X331VQYM\nGNA08bTStv6+Db+dhpYfpixHw5TJvP02DB8e/fzmmzBkSH3n6+joSPVLLOsot+Eot+Eot+Eot0Gp\ntYU0zs03w5o1UV+xegsx0FppISm34Si34Si34Si3+dCuxZh6xydQaGmhtShFRETSl2iY0sz6AJOB\nA4HdgKFAH2AxMAd4ALjB3eekHWgazGwWMIZozlhv4E2iOWV7uvuCkmM1TAm8+y5sthmsWAGvvALx\nlJm66LJ5OMptOMptOMptOMptUI29m9LMegL/CnwT2LRCIHsDJwLfN7ObgH9z95fSCjQN7r5f1jHk\nzZ//HBVie++dTiEmIiIi6+t2mNLMBgMzgYuBwaxfgHnRo/icxwKPmtkx6YUqWQjR6FW/pYWj3Iaj\n3Iaj3Iaj3OZDl8OUZtYPuBP4ELVdjlsNHO3uf64tvOxomBJWr4bNN4e33oKnnoJddsk6IhERkabR\nsGHK/wQ+HP/8HHAL8BgwF1gEvAesAfoCmwAjgZ2I5pRNJJqXda2Z7VI6J0ua39/+FhVi228PO++c\n3nk1hyEc5TYc5TYc5TYc5TYfKhZjZjYeOB14HviKu89McN6LzGxr4HvAJ4AfAJ+qJ1BpvOIhStP9\npyIiIkFUHKY0s5uJrnQd6u7v1PwBZpcQTf7fqdkm9Hel3Ycp3WHMmOgOylmzYN99s45IRESkqaR2\nmaJsMWbR2kGTgIfcfd4GByT9ELOjgKebteVFOe1ejD36KOy1VzRn7PXXoUe7dqQTEREpL2wHfo9M\nr7UQM7M+Zja66Hwz8lSIybpGr5Mnp1+IFa+VJulSbsNRbsNRbsNRbvMh1PWOfYDZgc4tDRCipYWI\niIhsqKaFws1sT6K7LIdQ/iaArYDPdrUQd7Nr52HKl16C7baDQYNg4ULo0yfriERERJpOYzvwr/1U\nswHADcBhaQUgzadwVeyoo1SIiYiIhJZ0mPL7wOF0Xw2qEUKOhR6i1ByGcJTbcJTbcJTbcJTbfEh0\nZQw4jmjpI0MFV0tasADuvTe6InbkkVlHIyIi0vqSFmNDgReBy4m68K8AOsscNzY+RnLmppuiHmOH\nHRbNGQtB3aDDUW7DUW7DUW7DUW7zIWkxNgf4mrvf1NVBZvY6sFetQUl2Ci0tjjsu2zhERETaRdI5\nY78CtqjiuMVEC4xLjixdCjNnRksfHXtsuM/RHIZwlNtwlNtwlNtwlNt8SFqMXQzsZ2ZfMLPN4079\n5ewAXFVfaNJot9wCK1fCfvvBiBFZRyMiItIeEg1TununmV0NTAeuBKhcj0neNKrRq+YwhKPchqPc\nhqPchqPc5kOipq/xGpPTiK6odVeFuZq+5sfKlTBsGCxZAi+8EDV9FRERkYrCrk3ZhX8HCgVWpUql\nfSqYFnLnnVEhtttu4QsxzWEIR7kNR7kNR7kNR7nNh6R3U+4WPy8EXqFya4tNgA/WEZc0mNaiFBER\nyUbSYcqlwI/c/RvdHDcBuMfdQy1EHlw7DVN2dsKoUTB/Pjz8MOylpiQiIiLdyWyY8kngr1Uc9wCw\ndfJwJAsPPBAVYqNHw557Zh2NiIhIe0lajP0E+GIXLS0Kdkd9xnKjuNFrI26O1RyGcJTbcJTbcJTb\ncJTbfEg6Z2wbYF9gtpndCiwA1pQ5bqv4IU3OfV0xpvliIiIijZd0zlgn6xYK7/aNam3R/J55BnbZ\nBYYMiRYJ75W0PBcREWlPmc0ZK+iuSlEn2JwoXBWbNEmFmIiISBZqKcasiofkRBYtLTSHIRzlNhzl\nNhzlNhzlNh9quRZyBtFdlV0ZC1xew7mlgV57DR58EPr3h8MPzzoaERGR9pS0GJsL3Ofuj3R1kJkt\nI2oKK01s2rTo+cgjYaONGve5WistHOU2HOU2HOU2HOU2H5IuFF5V7zB3fxj1GWt6xS0tREREJBtl\n54xZZKyZbZzGh5jZTma2WRrnknR0dMBdd0HPnnDMMY3+bM1hCEW5DUe5DUe5DUe5zYeyxVjc0+En\nwI1m1rueDzCzM4CngFH1nEfS9ac/wZo1cNBBUVsLERERyUZXd1NeCBwGzDSz7ZOe2MwGmtn3gP8E\n/uzu/6gxRgkgy4XBNYchHOU2HOU2HOU2HOU2H7ps+mpmfwCOI+qyfztwK/APoon8i4D33L3TzPoA\nmwAjgZ2BA4CTgKHASmCsu78Q8M+RulZu+rpiBWy2Gbz7LrzyCmy5ZdYRiYiI5E5qrby6m8D/RaJ1\nJrcFjowf60dSeTHDwgtfzlsh1ur+/OeoEBs3LptCrKOjQ7+tBaLchqPchqPchqPc5kOXTV/d/S3g\nI8DzrCuuqmn4akAncJa7X5V+2FKPLIcoRUREZH1VrU1pZoOAi4iulBXWmyz3xkIx9iRwhrvfnUaQ\nWWjVYcrVq2HzzeGtt+DJJ2HXXbOOSEREJJdSG6ZMulD4KOCTwIFEw5dDgD7AEmA28CBwA3BH3iuZ\nVi3G7r4bDj4Ytt8ennsOKo8yi4iISBcaNmdsPe7+GvD9+CE5VNzoNatCTHMYwlFuw1Fuw1Fuw1Fu\n86GWhcIlp9w1X0xERKTZJBqmzCMzGwCcAkwCdgE2Bt4F7gUucPcXK7yv5YYp//EP2HNPGDEC5s2D\nHirFRUREapXa+FI7/HM8DriSqDfaB919JNEdoh8EHjCzMdmF1liFIcrJk1WIiYiINIt2+Sd5nruf\n4e7LAeKrYecAmwKnZRpZAzXLEKXWSgtHuQ1HuQ1HuQ1Huc2HdijGHgGOKLP/1fh5kwbGkpmXX4bH\nH4eBA+GQQ7KORkRERAoS3U2ZR+6+FHi6zEvjiHql3dPYiLJRuCp29NHQt2+2sejOnnCU23CU23CU\n23CU23xo+WKslJltBBwOXAz8t7vfkHFIDVHc0kJERESaRzsMU65lZtcAi4HrgKnAv2QbUWO8+Sbc\ney/06QMf/WjW0WgOQ0jKbTjKbTjKbTjKbT6kWoyZWU8za9o5WO7+aaA/0QoCRwH/MLNtso0qvJtu\ninqMTZwIgwZlHY2IiIgUS1SMmdmBZjawi0N2AzrM7EEz27q+0MJw99Xufj9wArAN8IuMQwquMESZ\n9V2UBZrDEI5yG45yG45yG45ymw9J16ZcAxzg7rMqvD6YaDHxLwGPuPuJqUQZiJk9A2wPDHD390pe\n83POOYf+/fsDMG7cOCZMmLD2i1249Nvs2716DWbYMFi50nn66cXstNOmTRWftrWtbW1rW9s53c5s\nofBOYEKlYqzouMnAz919RJ3x1c3Mjgded/cHyrz2CFHz1y3c/Y2S11qiA//118NJJ8F++0XzxppB\nR4fWSgtFuQ1HuQ1HuQ1HuQ0qm4XCY11WKPHdinsAQ2qKKH3HAsuB9YoxMxsB7ATMLy3EWkmzNHoV\nERGR8rq8MmZm5wPfKt5FN8VYkRXuPqCO2FJhZlcRrU35f4H/dfdVZrYd8Etgf+DT7v7bMu/L/ZWx\nlSth+HBYvBheeAG22y7riERERFpGQ6+MlX5YtR/+csJYQjkPeIFoHtv3zKwvsAq4HzjU3Vu26etd\nd0WF2G67qRATERFpVtXeTVnLJaINrjZlwd1fc/fvufu+7r6Fu2/m7iPd/bhWLsSgeRu9FiY/SvqU\n23CU23CU23CU23zo7srYHODuou0DgceIGqeW0wksAm5z96l1Ryc16+yEadOinzVfTEREpHkFuZuy\nFeR9ztj998P48bDlljB3LlhqI9siIiJCinPGknbg/xzwfFofLuEUD1GqEBMREWleiYoxd7/a3Rd1\nd5yZbWdmGqbMUDO3tNAchnCU23CU23CU23CU23wItVD4COAzgc4t3XjmGXjuORgyBA44IOtoRERE\npCuJm76a2YeBKcC+wFCgZ9pBSX0KV8UmTYJetbT1DUzdoMNRbsNRbsNRbsNRbvMh0T/VZjYRmBG/\nr7uZSPmd/Z5zzdrSQkRERDaUdJjyAqB3N8eoCMvQa6/Bgw9C//5wxBFZR1Oe5jCEo9yGo9yGo9yG\no9zmQ9JBrD3j57eIepAtr3DcJkQLcEuDTZ8ePX/kI7DRRtnGIiIiIt1LWoytAX4M/Iu7d1Y6yMwm\nAC3d3b5ZFYYom/EuygLNYQhHuQ1HuQ1HuQ1Huc2HpMOULwB/6qoQiz0OHFJbSFKrjo5oPcqePeHo\no7OORkRERKqRtBi7CvhoFceNBs5PHo7U4+abYfVqOPBAGDo062gq0xyGcJTbcJTbcJTbcJTbfEg6\nTPkUcLqZ/Ra4BXgdWFXmuLHAQXXGJgk1c6NXERERKa+WtSmrfoO757YHWd7WplyxAjbbDN59N1qL\ncvTorCMSERFpaaktNlhPS1D1GWsiM2dGhdi4cSrERERE8qSW5ZCMFKtBSUeeGr1qDkM4ym04ym04\nym04ym0+1HJl7CrglW6OGQ18toZzSw1Wr17XX0zzxURERPKlljljE9x9VjfHTQDucfdQC5EHl6c5\nY/fcAwcdBNttB88/D6brliIiIqFlNmdsa2B+Fcc9EB8rDVA8RKlCTEREJF8SXbly97nuvrJ4n5kN\nNrNeJcetdPe5aQQoXXPPX0sLzWEIR7kNR7kNR7kNR7nNh5qGEc3sIDO7ycwWAwuJ+ophZreZ2clp\nBihde+wxmDMHRoyA8eOzjkZERESSSlyMmdnXgb8ARwEDWX/MdDzwP2Z2VTrhSXcKV8UmT4YeOZmh\np7XSwlFuw1Fuw1Fuw1Fu8yHRP99mtj9wIesKsNIZ7h8Avg580sxOqj886U6eWlqIiIjIhpJeSzmT\ndX3GNug35u7L3f0S4CfA51KJUCp6+WV4/HEYOBAOPTTraKqnOQzhKLfhKLfhKLfhKLf5kLQY2w94\njmix8GFAf8rf2nkbsGd9oUl3pk2Lno86Cvr2zTYWERERqU3S1hbDgDPc/TYAM6u09mR/YEg9gUn3\n8jpEqTkM4Si34Si34Si34Si3+ZD0ytgKYPuuDjCzfsBZ8bESyJtvwr33Qu/e0ZUxERERyaekV8Ze\nBL5nZnsBDwKL4/3HxpP7dwSOA0YCT6QWpWzgppugsxOOOAIGDco6mmQ6Ojr021ogym04ym04ym04\nym0+JC3G/gSMA06KHwXnFf1sRHdZ3lpfaNKVvDV6FRERkfKSrk05GHgKGFHYRVR4Fbe6MKAD2M3d\nq1k6qSk189qUS5fCsGGwciXMmwebb551RCIiIm0ntQUIky6H1AEcTbQ+ZXEQharFgLeBj+W5EGt2\nt90G778P++6rQkxERCTvEvdsd/dHgZ2AfyNqYfEsUbuLmcA3gR3c/e40g5T15fUuygL1vQlHuQ1H\nuQ1HuQ1Huc2HpHPGAHD3ZcBl8UMaaOVKuPnm6Oe8FmMiIiKyTpDVDM1sNzO7I8S5293dd8PixbDr\nrrB9l01Gmpfu7AlHuQ1HuQ1HuQ1Huc2Hmq6MAZhZf2BwhXPsCBxU67mlsrwPUYqIiMj6El8ZM7NT\nzOwxYCnwKjC7zOO6NIOUSGfnuiWQ8tzSQnMYwlFuw1Fuw1Fuw1Fu8yHRlTEz+xRwdWGzm8Obsy9E\njj34YNTKYsstYa+9so5GRERE0pB0mPJrpNhXQ5IpNHo97jiwHP9X0ByGcJTbcJTbcJTbcJTbfEha\njO0AdAIPA3OJ1p/sLHPccODI+kKTUpovJiIi0nqSduB/E7jI3btsaWFmE4B73D3I3ZqN0Gwd+J99\nFnbeGQYPjhYJ71XzrRfZ01pp4Si34Si34Si34Si3QaU2RpX0n/Q/Ey0W3iV3/xuB2ma0q8IQ5aRJ\n+S7EREREZH1JC6ZvAqeZ2cCuDjKz/c1sde1hpcfMBpnZmWZ2n5ktMrMOM3vCzL5mZrkpa1ppiFK/\npYWj3Iaj3Iaj3Iaj3OZDomFKiBq6An8BHgUWAGvKHDYc+Ki796w7wjqZ2QyinmenuPt0M+sJfAb4\nOXCzu0+u8L6mGaZ8/XUYNQr694dFi2CjjbKOSEREpO1lM0xpZpsT9RAbBhze1aE0T2sLAy539+kA\n7r4GmGpmE4GTzWyiu/8l0wi7UegtdsQRrVGIaQ5DOMptOMptOMptOMptPiQdpryEaJHwZim0qvE/\nwDVl9t8XP+/TwFhqUpgvludGryIiIlJe0rsp3yAagqyGN8MwZSVmdjZwKfBP7n5lmdebYpiyowOG\nDwd3WLAAhg7NOiIREREhw7spNwX+DnyXrvuM7UPzL4m0D7AKmJ51IF2ZMQNWr4ZDDlEhJiIi0oqS\nFmPzgO+6+y1dHWRmncB3ao4qMDPbEjgWuMLd52cdT1eKu+63Cs1hCEe5DUe5DUe5DUe5zYekc8Zu\nBHav8tjsx/gq+xnwJPDvWQfSlRUr4Ja47G2lYkxERETWSXpl7DfAtLjP2M1Ubm2xD3A+cEF94aXP\nzC4luglhvLuv7OrYKVOm0L9/fwDGjRvHhAkT1v6G0dHRARB0+7bberN8+QD22gsGDuygoyPs5zVq\ne/DgwU0Vj7a1Xe12QbPE0yrbhX3NEk8rbevv28Z8f+uVdAJ/JwmueDXbBH4z+zrwz8AB7j6nm2Mz\nn8D/+c/D1KlwwQVw3nmZhiIiIiLrS20Cf61LFlkVj6ZiZmcA/wJMLBRiZjbEzLbKNLAK1qyB6fGt\nBa3W0qL0KoOkR7kNR7kNR7kNR7nNh7ZYP9LMTgPOAw539+eLXppENJzadO69N+q2v+22sOuuWUcj\nIiIiodSyNuOtwJvdHDMcOLKGc6fOzE4mWvroT8AJZnZC0ct7AE35a0Nxo1druuuM9UlznF3Wp9yG\no9yGo9yGo9zmQy1zxia4+6xujpsA3OPumV95M7NHgbFdHPIrd/98mfdlNmfMHbbZBubMgb/9Dfbf\nP5MwREREpLLM5owdAjxRxXGPx8dmzt33dPeeXTw2KMSy9vjjUSE2YgSMH591NOnTHIZwlNtwlNtw\nlNtwlNt8SDRM6e53V3ncEqCqY2VDhSHKY4+Fnk11P6qIiIikraZhRDPb3sz+y8yeMLN3zGzPeP/P\nzawFr+VykVl2AAAgAElEQVQ01o03Rs+t2uhVcxjCUW7DUW7DUW7DUW7zIdGcMQAz+xTRhPi+ROOl\nDuzj7o+Y2XtAb6Ilk76dcqwNldWcsdmzo/liAwZEd1P27dvwEERERKR72cwZM7PdgV8SFWLl7Eq0\n1NB5ZnZEnbG1pcIQ5VFHtW4hpjkM4Si34Si34Si34Si3+ZC0tcVZRFe+ynL3l4B/NrONgS8Bt9cR\nW1sqbmkhIiIirS9pa4uXgY2BbwEPAcuAZ4C93f2RouMmAT9x9y3TDbdxshimXLgQNt88mrS/cCFs\nsklDP15ERESql9owZdIrYyOBz7r77wDMrNK9fp1EjV8lgZtugs5OOOIIFWIiIiLtIundlKuAgVUc\nd3p8rCRQGKJs1bsoCzSHIRzlNhzlNhzlNhzlNh+SXhmbA1xuZtsCDwJL4v17m9kWwI7AycCewHNp\nBdkOli2D22+Plj6aPDnraERERKRRks4ZuxQ4m6idxdrdFbZ/6u5fSSPILDR6ztgNN8CJJ8K++8Ks\nLhebEhERkSaQ2XJI/8m6q2HGusLLih4A7wNXpBFgu2j1Rq8iIiJSXqJizN1fB04iuosS1l0R86Kf\n3yea5P9CKhG2gVWr4E9/in5uh5YWmsMQjnIbjnIbjnIbjnKbD4mXQ3L324HdgR8BzwPvERVgLxF1\n5t/D3a9LM8hWd9ddsHgx7LILbL991tGIiIhIIyWdwA+Au78CnJlyLG2r3Rq9aq20cJTbcJTbcJTb\ncJTbfKhpofBKuug7JhV0drZPSwsRERHZUNK1Kbc2szuKHhuVHPKwmf3czDZNMcaW9tBDMG8ejBoF\n48ZlHU1jaA5DOMptOMptOMptOMptPiS9MnYUcDBRP7Fn2LCx6wvAqcAsFWTVKb4qZqndJCsiIiJ5\nkbTP2DRgc2Ciuy+rcMx2RAuEX+3u30klygw0qs/YzjvDs8/CzJkwcWLwjxMREZF0ZNZnbHvgskqF\nGIC7vwhcCLTJdPTaPfts9Bg8GA48MOtoREREJAtJi7HNgYVVHDcb2Cp5OO2lMER5zDHQu3e2sTSS\n5jCEo9yGo9yGo9yGo9zmQ9JirD9wUBXHHQD0SR5Oe2m3lhYiIiKyoaRzxmYTXR07D/hfd59X8voo\n4FPAt4HX3H279EJtrNBzxl5/PbqDsl8/WLQINt442EeJiIhI+lKbM5a06ess4BTg+8D3zWwZsDh+\nbVNg45JjpYLp06Pnj3xEhZiIiEg7SzpMeWXRzwYMBEbFjwGsv1j4lUhF7dzoVXMYwlFuw1Fuw1Fu\nw1Fu8yHpQuH3AJcTFVxe4QHwQ3e/N8U4W8o778Add0CPHtHkfREREWlftSwUfjbwr8BbrLsSVni8\nDZzt7melGWSrmTEDVq+O2llstlnW0TSe1koLR7kNR7kNR7kNR7nNh1oXCv+hmf0Y2BsYHe9+BXjI\n3dekFVyruvHG6LkdhyhFRERkfYnupqz6pGaDgD3iYc1cCnU35XvvRVfDli+HOXNgqzbsxtbR0aHf\n1gJRbsNRbsNRbsNRboPKrAN/tXYH7gx07lybOTMqxPbcsz0LMREREVlfTcOUZrYjcBxRl/1+ZQ4Z\nXk9QrUyNXjWHISTlNhzlNhzlNhzlNh8SD1Oa2TeAC+j6qpoB7u4964gtUyGGKdesgZEjYeFCePxx\n2H33VE8vIiIijZPNMKWZHUC0CHjhfakF0g5mzYoKsW23hd12yzqa7KjvTTjKbTjKbTjKbTjKbT4k\nHaY8PX4uFGGOCrKqFTd6NWVNRERESL425fPAYOCfgceAZ4BJwJPxIYOAQ4iKto+7+zOpRttAaQ9T\nukdXxGbPhr/+FSZMSO3UIiIi0niZrU05ErjA3a8zs8J8sIXuPrfomCfMbCXwMaIhTQGeeCIqxIYP\nh333zToaERERaRZJW1v0Af4BEDd37SRq/FpqAfCJ+kJrLYVGr8ceCz1ze1tDOjSHIRzlNhzlNhzl\nNhzlNh+SXhlbVPKeDuB7ZtYDuBtYBewKXAR8IJUIW4RaWoiIiEg5SeeMPQzc6+5nxtt3AAezboHw\ntYcCb7r75inF2XBpzhmbPRu22QYGDIjupuxXrjObiIiI5ElmHfhnAf9sZpfF29cWvVa8YLgDf68/\nvNYwbVr0fNRRKsRERERkfUmLsT8A9wDD4u1fArezYXX4JvDV+kJrHcUtLURzGEJSbsNRbsNRbsNR\nbvMh0Zwxd7+TojUn3X2NmX0UOBGYQHRF7DngN+6+NM1A02BmnwZ+CNzo7qc14jMXLoxaWfTuHV0Z\nExERESlW09qUxeKJVdfHj6ZkZkOBnxHd+TmokZ/9pz9BZyccfjhsskkjP7l5aa20cJTbcJTbcJTb\ncJTbfEg6TFkVM9vfzNaEOHeNfg28ABxBg1cMKLS00BCliIiIlFPxypiZnVrHeXeo470hfNHd55nZ\nVo380GXL4Pbbo58nT27kJze3jo4O/bYWiHIbjnIbjnIbjnKbD10NU/6KDVtW5JK7z8vic2+7Dd5/\nH8aPh5Ejs4hAREREml01c8aSDusVFg9viUKuHmr0Wp5+SwtHuQ1HuQ1HuQ1Huc2H7uaM1TK/qqFz\nsprVqlXR5H3QfDERERGprLsrYxcCL9Zw3h2Ar9fwvpZx993wzjuwyy6wQ7PNoMuY5jCEo9yGo9yG\no9yGo9zmQ3fF2C3uPivpSc1sAvCN2kJqDWr0KiIiItWoWIy5e81tL9z9bwRqm9FIU6ZMoX///gCM\nGzeOCRMmrP0No9DVuNx2Zyf84Q+dQI+1xVhXx7fb9uDBg5sqHm1ru9rtgmaJp1W2C/uaJZ5W2tbf\nt435/tYr0ULh673RrC8wDtgi3vU68LC7v59SbKmLW1vMBn7VXQf+ehYKf/BB+NCHYIst4NVXwTSL\nTkREpNVktlA4ZtbXzC4GFgJ/JVos/Nr45zfN7Ptxoda2ihu9qhDbUOlVBkmPchuOchuOchuOcpsP\niZZDiousPwP7U74iHAj8GzDezA5395X1h5iqhpRGamkhIiIi1Uo0TGlm5wAX0XUPsUKPsa+7+yX1\nhZcOMzsFuIzoSuBmwHvAEmCBu+9R4T01DVM+9xzstBMMHgwLFkQLhIuIiEjLyWyY8v8UBWDAO0Rz\nxebFP1vR659MI8A0uPtv3X2ku49w957uvnG8XbYQq0fhqtgxx6gQExERke4lLca2B+YAxwID3X2o\nu4929y3dfSjRMOVkYG58bNtRS4vuaQ5DOMptOMptOMptOMptPiSaMwZ0Aue4+5/Kvejuy4GbzKwf\nMLXe4PJm3jz4+9+hXz/4yEeyjkZERETyIOmVsaeBnlWe9+niHWa2nZm1dIE2fXr0fMQRsPHG2cbS\nzNLszSLrU27DUW7DUW7DUW7zIWkx9gPgfDPbqdIB8WvfAi4teWkE8JmEn5crxS0tRERERKqR9G7K\nO4kavW5MNHH/VaDQ5LUfUQPYUcAy4OGSt28C7OHu1VxZy1zSuynfeQeGDYPOzuguys02CxhczhV3\n2pZ0KbfhKLfhKLfhKLdBpXY3ZdI5YwcRta0woqJri5LXC4ENiI8tfa22lvY5MGMGrF4NBx2kQkxE\nRESql7QYK6hUVLVssdUdNXqtnn5LC0e5DUe5DUe5DUe5zYekw5SddX6et+Iw5XvvRUOUy5bB7Nkw\nZkzY2ERERCRz2a1NCUxw9x5JH8CBaQXdbP7yl6gQ23NPFWLVUN+bcJTbcJTbcJTbcJTbfEhajF0N\nLKjxs94Afl3je5uaGr2KiIhIrRINU7aTaocp16yBkSNh4UJ47DEYO7YBwYmIiEjWshmmNLOtzez0\n+DE+3tfTzL5lZi+b2XIzu8/MWnZIstR990WF2DbbwO67Zx2NiIiI5E3SYcrPAD8DfgoUFvz5FvBt\nYAzQH/gwcIuZtcXalMWNXi21Grm1aQ5DOMptOMptOMptOMptPiQtxvYmavY61t2/Y2YbA2cRtbQo\nPCBqAPuV1KJsUu5qaSEiIiL1Sdra4jngJ+5+Rbx9AnA964qwhcAiYBfgcXffI91wG6eaOWOPPw4f\n/GDU1mL+fOiZi6YdIiIikoLMWlsMA54s2j48fjbgTWAHd98N+A3RsGVLK1wVmzxZhZiIiIjUJmkx\nNgDYqGh7YvzswHXuviTe/l3JcS1JLS1qozkM4Si34Si34Si34Si3+ZC0GFsGnAJgZkcA2xa99uei\nn1cDK+oLrbnNmQOPPgoDBsDEid0eLiIiIlJW0rUpXwA+YWYfAQaxbtHw94E7i44bB7ydSoRNatq0\n6PmjH4V+/bKNJW+0Vlo4ym04ym04ym04ym0+JL0yNoOo+NoUKMyScmCGuy8HMLNxwNnAy2kF2YyK\nW1qIiIiI1CppMXY58CLr7iAwYClwDoCZ3Qj8HRgCPJxSjE1n0SL461+hd284+uiso8kfzWEIR7kN\nR7kNR7kNR7nNh0TDlO6+OL7y9XlgR6KeY1e5++vxIfcAz8Y/35BalE3mppugsxMOOww22STraERE\nRCTPtDZlBV31GZs8GaZPh5/+FL70pQYHJiIiIs0gtT5jQYqx+OrZ9e6+Teonb5BKxdjy5bDZZvDe\ne/D66/CBD2QQnIiIiGQttWKsy2FKM9sGwN1fjrerXQB8LLBVfaE1p9tuiwqx8eNViNWqo6NDd/gE\notyGo9yGo9yGo9zmQ8VizMy+AfxH/PM33f37wF2sW/qoLanRq4iIiKSp4jClmXUAhenpHe4+1Mw6\nWddbrDvu7rldJKjcMOWqVTB8OLzzDjz7LOy4Y0bBiYiISNYasjbl/Ao/p/bheXPPPVEhtvPOKsRE\nREQkHV3NGfs0cAXRlbB/Kdp/FfBKN+cdDXy2rsiakBq9pkNzGMJRbsNRbsNRbsNRbvOhYjHm7g8D\nE8q89Et3n9XVSc1sAvC5OmNrKu7r5osdf3y2sYiIiEjrSNTawsy2Aua7+8pujusDjHT3uXXGl5nS\nOWMPPQT77ANbbAGvvAI9kq5dICIiIq2kMa0tShWKKzMbDQyLdy9091dKjlsJ5LYQK6d4iFKFmIiI\niKSl6rLCzMaY2S/M7A1gNvBA/JhtZvPN7KdxkdaS1NIiPVorLRzlNhzlNhzlNhzlNh+qKsbM7OPA\n08BpwHCiS3PFjxHA6cAzZtZyM6qefx6efho23RQOOijraERERKSVdFuMmdmhwLVAP6LCyys8DOgP\nXJegU38uFK6KHXMM9O6dbSytQHf2hKPchqPchqPchqPc5kOXxZiZ9QSmsq4Igw2vihnrF2k9galm\n1jIzq9TSQkRERELp8m5KMzsBuJ51V77+AdwHvAQsifcNBLYF9gM+GL/VgY+7+x+CRR5Y4W7K+fOj\nNSj79YNFi2DjjbOOLP/U9yYc5TYc5TYc5TYc5Taoht1NeXT8/CTwOXd/pKuDzWwcUVPYXYFjgNwW\nYwXTpkXPhx+uQkxERETS192VsSeAzYBd3L2qWzLMbCjwFPCmu49NJcoMFK6MHXkk3HYbTJ0Kn2up\nNrYiIiJSh4asTQmwOXBltYUYgLu/BVwJjKwnsGaweDHccUfUV+yYY7KORkRERFpRd8XYpkS9xJJ6\nANikhvcFY2YnmtlDZrbAzOaa2SVm1r+r98yYAatWwYQJMGxYV0dKEup7E45yG45yG45yG45ymw/d\nFWM9gXdqOO878XubgpmdBvwOuNTdRwAHApOBm8ys4mVGrUUpIiIioXU3Z6wT+ArRBP4kxgJXuHvm\nBZmZbUq0YsCt7n5K0f5JwDSiGxOuLvM+HzDAWbYMZs+GMWMaFrKIiIg0v9TmjFVTjFW/kniJJinG\n/i/wE+AUd7+uaH8vYDHwiLsfUOZ9Ds4ee8CjjzYuXhEREcmFhk3gL/7ApI9mUSi0Hi/e6e6riZZ4\nGm9mFfvqq9Fr+jSHIRzlNhzlNhzlNhzlNh+qLcYqLYFU6dFMdoif55d5bR5RDrap9GbNFxMREZGQ\numv6CrAAeD/hefsSLR7eDAp3db5b5rXCvk3LvXHrrWH33YPE1NbUDToc5TYc5TYc5TYc5TYfqinG\nTnT3WUlOamYTgLtrC6l5HH88VL7XUkRERKR+3RVji4HVNZx3FdHalc1gcfy8UdHPFO2jzH4ALhtk\nXPad9ff5+eVHYe075as2Ha/jG3n822e+XfY34bzE38zHF6/x1wzx6Hgdr+Ob8/hadDlnzN0Hu3vi\npq/ufr+7N8u10efj53IrAnwA6AReTnrSjo6ORBMjdbyO1/H5Pn7ZsmVNFY+O1/E6vjmOT0OXrS1a\ngZmdDvwMOLlMa4t3iFpbHFjmfd7quREREZGaNby1RZ5dTzRkWnpf5FFEw5S/bHhEIiIiIrGWL8bi\nRc7PBk40s1MAzGwMcAnwF+DXmQXXptT3JhzlNhzlNhzlNhzlNh9avhgDcPepwMnAv5nZAuAeoqWQ\nJmksUkRERLLU8nPGaqU5YyIiItIFzRkTERERaQUqxqThNIchHOU2HOU2HOU2HOU2H1SMiYiIiGRI\nc8Yq0JwxERER6YLmjImIiIi0AhVj0nCawxCOchuOchuOchuOcpsPKsZEREREMqQ5YxVozpiIiIh0\nQXPGRERERFqBijFpOM1hCEe5DUe5DUe5DUe5zQcVYyIiIiIZ0pyxCjRnTERERLqgOWMiIiIirUDF\nmDSc5jCEo9yGo9yGo9yGo9zmg4oxERERkQxpzlgFmjMmIiIiXdCcMREREZFWoGJMGk5zGMJRbsNR\nbsNRbsNRbvNBxZiIiIhIhjRnrALNGRMREZEuaM6YiIiISCtQMSYNpzkM4Si34Si34Si34Si3+aBi\nTERERCRDmjNWgeaMiYiISBc0Z0xERESkFagYk4bTHIZwlNtwlNtwlNtwlNt8UDEmIiIikiHNGatA\nc8ZERESkC5ozJiIiItIKVIxJw2kOQzjKbTjKbTjKbTjKbT6oGBMRERHJkOaMVaA5YyIiItIFzRkT\nERERaQUqxqThNIchHOU2HOU2HOU2HOU2H1SMiYiIiGRIc8Yq0JwxERER6YLmjImIiIi0AhVj0nCa\nwxCOchuOchuOchuOcpsPKsZEREREMqQ5YxVozpiIiIh0QXPGRERERFpB2xRjZnakmb1mZndkHUu7\n0xyGcJTbcJTbcJTbcJTbfGj5YszMNjKzHwM/BUZmHY+IiIhIsZafM2ZmPwN6Af8KLAHucvdDq3if\n5oyJiIhIJanNGeuV1oma2AXuPg/ALLW8iYiIiKSi5YcpC4WYNA/NYQhHuQ1HuQ1HuQ1Huc2Hli/G\npPn87W9/yzqElqXchqPchqPchqPchmNmB6d1LhVj0nAPP/xw1iG0LOU2HOU2HOU2HOU2qIPTOlEu\n5oyZ2UbA8QnecqO7vxsqHhEREZG05KIYA4YB1wDV3t64PfByuHCkHitWrMg6hJal3Iaj3Iaj3Iaj\n3OZDblpbmNmgao919yUVztFJgtYWCcITERGRNuPuqbRpyMuVsYoFVsDPUx8MERERCU4T+EVEREQy\npGJMREREJEMtX4yZ2QQzm29m84huANivsG1mvUuOPdHMHjKzBWY218wuMbP+2UTeWsxsTpzz4sd8\nM3sl69jyxsw+bWYdZja1i2OGmtkv4zy/YWZ3m9lBjYwzj7rLrZl9xsyWlvkezzOzMxsdb7Mzs0Fm\ndqaZ3Wdmi+LcPmFmXzOzDabJ6HtbvSS51fc2GTMbYGZfNLPpZvZinKuXzOzXZrZdmePr/966ux7R\nTQynAWuAk+PtrYDngZnENzroUVd+X846hrw/gKHA9cDs+Ls6tcJxA4CngL8Cg+N9XwdWAROz/nM0\n4yNBbj8DfCvrePPyAGYAy4Fj4+2e8d+1q4FpJcfqexsut/reJsvtQUAn8F/AxvG+7YDHgLeBMUXH\npvK9bfkrY9Uws02BHwDXufu1AO4+F/gqcChwaobhiRT8GngBOIKuF6g9B9gJ+IK7dwC4+0VEf5H8\nzMz0//2Gqs2tJGPA5e4+HcDd17j7VOB3wDFmNrHoWH1vk0mSW0lunruf4e7LAdz9RaLv6KZERW9B\nKt9bfbkjnwAGATeW7L8FWAF8oeERiWzoi+7+TWBlN8edBjzn7s+V7P8DsA1wSIjgcq7a3Eoy/0PU\nI7LUffHzPkX79L1NJkluJZlHiH4xK/Vq/LxJ0b5UvrcqxiIHxM+PF+9099XA08D40vllIo3mVSx6\nH89n+AAl3+XYY0S/TWsOTolqcivJuftv3P3ZMi/1IfoudoC+t7WoNreSnLsvdfeny7w0jmju+T2Q\n7vdWxVhkh/h5fpnX5hHlaZvGhdOSzMwuNLMn48mQT5vZD8xsaNaBtZjuvssQrVAhtdvHzG4puinl\ndjOblHVQObMP0Zya6fG2vrfpKc3t2v363tbGzDYys8nAxcB/u/sN8UupfW9VjEUKlxzLrWdZ2Ldp\ng2JpVZ1EudwXGAV8Gfg48KCZDc8ysBaj73J4Y4Dz3H0M8EHgOWCamZ2TZVB5YWZbAscCV7h74R8x\nfW9TUCG3BWPQ9zYxM7sGWAxcB0wF/qXo5dS+tyrGpFH2cfcL48u/a9z9LqKCbAzwH5lGJlK93xF9\nlx8CcPeF7n4G8DDwXTMbnWl0+fAz4Eng37MOpAVVyq2+tzVy908D/YEDgaOAf5hZ6iNlKsYii+Pn\njcq8tlHJMVIDd3+7zO4ZRLdhH9PgcFqZvssBuft77v5emZduIlpe7sgGh5QrZnYp0Z1nx7h78c0S\n+t7WqYvc6ntbJ3df7e73AycQTVn6RfxSat9bFWOR5+PnkWVe+wDRENvLjQunPbh7J/AWMCzrWFpI\nd99liFo4SLoWEE3W1ZB7BWb2daI71ye6+8KSl/W9rUM3ue2KvrcJuPvLwEvAQWbWjxS/tyrGIvcQ\nfSHHFu+MuxjvDNxX+puGVM/MDjKzw8rs70HUbPOtxkfVmuJeOPMo+S7HxhLdCXRXI2NqJWZ2vpn1\nLPPSCKLcLmpwSLlgZmcQzbWZ6O5z4n1DzGwr0Pe2Ht3lNt7W9zYBMzvezD5U4eUVRPXCpml+b1WM\nRa4HlgDHl+w/iuhS4y8bHlFrORg4o8z+I4kukd/S0Gha31RgRzPbqWT/iUS/1d3Z+JBaxvnA7mX2\nH010Bf32xobT/MzsNOA84HB3f77opUlE+SzQ9zahbnL77aJtfW+TOZYyzd7NbATRUPB8d38j3p3K\n91bFGBB3zT0bONHMTgEwszHAJcBfiLpzS32OMbMvF/q1mdm+wI+Ibgk+L9PI8qe7DvEXA88AP4/X\nTDMz+wbRX8b/FA8PS3nVdN//mZntAGvXB7wU2Bv4fjyMITEzOxn4OTALOCG+QnO+mZ0PHFdyuL63\nCVSRWy95i763yXwhXtOz8G/WdkR3VPYBvlZ0XCrfW4vXURLAzE4AvknUeuF94Fqi9bzKTXyUKsW9\nxD7FusmP/YGlRBP4v1vmFmwpI/5F4TKiX6I2A94juqK7wN33KDl2CPB9ot98exDNbfh3d7+noUHn\nRLW5NbODgU8D+wGDgX7AE8CPC0upyTpm9ijlh3AKfuXuny86Xt/bKiXJrb63yZjZKKIrY5OA0UBf\not5t9wOXlX4f0/jeqhgTERERyZCGKUVEREQypGJMREREJEMqxkREREQypGJMREREJEMqxkREREQy\npGJMREREJEMqxkREREQypGJMRKTJmFk1KwGk/l4RyYaKMRFpevHCx382s6Vm9rvCEiWtJv5zXky0\n3l2tvmdmPzWzLdOKS0TCUjEmIsGY2V1m1lnHY3R8qrOAicBGRAvwfjqzP1QgZvYJ4Dngy8D0Ok41\ngyhHT5vZV9KITUTCUjEmIiE58CiwPzCh6PFG/Frh8ZWS17/HhgsdQ3ULeeeOmX0L+C0wBPi4u99Y\n67nc/a/A4UAv4IdmdmU6UYpIKL2yDkBEWt5id7+veIeZvV/4kajoetLdZxW9vn3JOS4DPgyMJ7ry\nc024cBvLzL4EfJsoD//j7rfUe053/4eZ/QD4JvAFM1vk7ufWe14RCUNXxkQkpNnxI6kO4ClgJYC7\nd7j7Ee4+yN1PdvdVaQaZlXgY9jLWXQX8rxRP/6Oi804xs71TPLeIpEjFmIgE4+6nufvna3jfdHcf\n6+5vmNnsMnPJvlU41syuKvP6HWa2r5ndZGZvmdkKM3vGzM4uet/HzOw+M1tmZh1m9kcz26FSTGa2\nt5ldbWYvx+dbYmaPm9n/M7NhybMDwL8D/YiuEL7t7g9W+OzDzOy6OBfvxp//mpk9YGb/bWafLJPD\nN4DH4nMbcEGNMYpIYOZeblqGiEg4ZjYbGM26YcpD3P2eCsfuBexLdNWo8BfWd9z9gvj1bYERwN+K\nXu8AXgO+C2wSv7d//Pp3gKXAcOBuognzx8TvewPY1d07SmL4BvAfcbwvAFOA7YCLiH6pXQAc4e5P\nJMhBn/h9g+Jdd7n7xDLH/RPw43hzQRz/c8CmwMeBk4E33P0DZd57FfCZeHMNMMrdF1Qbo4g0huaM\niUhTc/dHzGxgF6+/BLxU1F7LiAqV8e7+IoCZfYSocAH4BnCdu38mfm0uUTHmREXdZ4H/XHsys48B\nF7LuZoMz3P32+LVxwCeICrvfm9ku7r6myj/aB4kKxUIBObf0ADPrQVQEFo75f+7+s6JDbjSzXkQ3\nPZRTfM4e8XE3VBmfiDSIhilFpBW9WCjEYoWixIA+wC1lXiv4cMn2+UTFUKHae6joteKftwOOSxDj\n7iXbi8ocMwwYXLR9eJnC9AfAzyt8RuGchWJutwTxiUiD6MqYiLQaJxrOK7aiZHtBF68NKfxgZpsT\nFU3F8zmWmVnP+OelhUPjYw6m+itPQ0q23ytzzEKimxgKTW6PBhaY2T3AHcDN7n4/cH+Fzyg959Aq\nYxORBtKVMRFpRaVDhZ2VXvcNJ84W/5Ja2sXeiAqcVfHjp6wbvoRoHly1+hSds1yMuHsn0V2RxfoC\nRxDNV3sivglhzwqfUZqHPmWPEpFM6cqYiEhlpU1mHfgQ5RvSAixJcO7CsYVzlS2U3P3f4hsezgC2\nL/l/UosAAAKmSURBVHpP4Wrch4FbzWxnd3+75O19S7YXJ4hPRBpExZiISGWvldn3vLsvLbM/qdL+\na5tUOtDdfwz82Mx2AQ4FJsfPBZsR3YTw6wrnLBRuc+qIV0QC0TCliEgF7j4PeJ71r5DtWHqcmZ0W\n9/76VILTP1SyPbLMeXvH/dEOjeN52t1/5O6HA+eyrsiC6E7QUqXnfCBBfCLSICrGRES6dlHJ9noN\nVs2sN9HamkuBqteUjPt93cu6pqw7lzmsB1Hxd2aZ1wp3ixYKxXI9znYq+nmuuz9cbXwi0jgaphSR\nhjGz/YmKh76sf7VprJmtIRoCXFjynnHA2MIm0ZWgreJzPUg0aX7zotcBNil6fSwbTqwfa2bziOZQ\n7Vjhvc+7+0J3/5WZ7Q78a/z6mWa2GriVqPXEV4iKnuPdfXmyjPBDokXUAbY3s2Glf/7YJDObAVxN\ndIflDkTrThauis1w91uL3xD3H9u36JgrEsYmIg2iDvwi0jBm1knlye8An3P39eY9FXXrL2droo70\np3bx+t0V3n91/NpVFWJaLxYzOwT4ElGBMyx+z2vAXcBl7v5shRi6ZGYzieZ/OXCeu3+v6DUjKgIL\nBelwooa2TtRD7AngeuDq+M7L4vN+inVzyB4HxiVoSCsiDaRiTEQkQ2Y2FLiPqGnsEmA7dy/XADbJ\nOfsAzwJjiJZ42tfdN+jwLyLNQXPGREQy5O5vEQ1VzgIGAn8xsw0m81cr7tB/C7AV8BQqxESanoox\nEZGMxfPEDgDOAragaG3MGpwH7ANcQDQ0qUJMpMlpmFJEpImYWT+iIureGt+/N/CCu6vBq0hOqBgT\nERERyZCGKUVEREQypGJMREREJEMqxkREREQypGJMREREJEMqxkREREQypGJMREREJEMqxkREREQy\n9P8B5Uzn1XCEWGcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113027080>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set the plot size - 3x2 aspect ratio is best\n",
    "fig = plt.figure(figsize=(6, 4))\n",
    "ax = plt.gca()\n",
    "plt.subplots_adjust(bottom=0.17, left=0.17, top=0.96, right=0.96)\n",
    "\n",
    "# Change the axis units to serif\n",
    "plt.setp(ax.get_ymajorticklabels(), family='serif', fontsize=18)\n",
    "plt.setp(ax.get_xmajorticklabels(), family='serif', fontsize=18)\n",
    "\n",
    "# Remove top and right axes border\n",
    "ax.spines['right'].set_color('none')\n",
    "ax.spines['top'].set_color('none')\n",
    "\n",
    "# Only show axes ticks on the bottom and left axes\n",
    "ax.xaxis.set_ticks_position('bottom')\n",
    "ax.yaxis.set_ticks_position('left')\n",
    "\n",
    "# Turn on the plot grid and set appropriate linestyle and color\n",
    "ax.grid(True,linestyle=':', color='0.75')\n",
    "ax.set_axisbelow(True)\n",
    "\n",
    "# Define the X and Y axis labels\n",
    "plt.xlabel('Time (s)', family='serif', fontsize=22, weight='bold', labelpad=5)\n",
    "plt.ylabel('Displacement (m)', family='serif', fontsize=22, weight='bold', labelpad=10)\n",
    "\n",
    "# Plot the data\n",
    "plt.plot(sim_time, response[:, 0], linewidth=2, linestyle='-', label = r'$x$')\n",
    "plt.plot(sim_time, response[:, 1], linewidth=2, linestyle='--', label = r'$\\theta$')\n",
    "\n",
    "# uncomment below and set limits if needed\n",
    "# plt.xlim(0, 5)\n",
    "plt.ylim(-1, 5)\n",
    "\n",
    "# Create the legend, then fix the fontsize\n",
    "leg = plt.legend(loc='upper right', ncol = 2, fancybox=True)\n",
    "ltext  = leg.get_texts()\n",
    "plt.setp(ltext, family='serif', fontsize=20)\n",
    "\n",
    "# Adjust the page layout filling the page using the new tight_layout command\n",
    "plt.tight_layout(pad=0.5)\n",
    "\n",
    "# Uncomment to save the figure as a high-res pdf in the current folder\n",
    "# It's saved at the original 6x4 size\n",
    "# plt.savefig('SpringBar_Pendulum_InitCondResponse_CartAndRadial.pdf')\n",
    "\n",
    "fig.set_size_inches(9, 6) # Resize the figure for better display in the notebook"
   ]
  },
  {
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}