CDPM_ODE.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint\n",
    "import numpy as np\n",
    "from numpy import sin, cos, sqrt\n",
    "import seaborn as sns\n",
    "sns.set_context(\"notebook\", font_scale=1.5, rc={\"lines.linewidth\": 2.5})\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def eq_of_motion(w, t, p):\n",
    "    \"\"\"\n",
    "    Defines the differential equations for the coupled spring-mass system.\n",
    "\n",
    "    Arguments:\n",
    "        w :  vector of the state variables:\n",
    "                  w = [x, theta, x_dot, theta_dot]\n",
    "        t :  time\n",
    "        p :  vector of the parameters:\n",
    "                  p = [m1, m2, k, l, g, wf]\n",
    "    \n",
    "    Returns:\n",
    "        sysODE : An list representing the system of equations of motion as 1st order ODEs\n",
    "    \"\"\"\n",
    "    x, y, beta, x_dot, y_dot, beta_dot = w\n",
    "    m, k, g, H, c, D, t, Izz, k_beta, c_beta, L_1_init, L_2_init = p\n",
    "\n",
    "    # Create sysODE = (x', theta', x_dot', theta_dot'):\n",
    "    sysODE = [x_dot,\n",
    "              y_dot,\n",
    "              beta_dot,\n",
    "              (-D*m*sin(beta)*beta_dot**2/2 + D*m*(-D*g*m*sin(beta)/2 + D*(-D*m*sin(beta)*beta_dot**2/2 + c*(H - x)*(x_dot + y_dot)/sqrt((H - x)**2 + y**2) - c*(x_dot + y_dot)*x/sqrt(x**2 + y**2) - 1.0*k*(-H + x)*(-L_2_init + sqrt((H - x)**2 + y**2))/sqrt((H - x)**2 + y**2) - 1.0*k*(-L_1_init + sqrt(x**2 + y**2))*x/sqrt(x**2 + y**2))*cos(beta)/2 + D*(D*m*cos(beta)*beta_dot**2/2 - c*(x_dot + y_dot)*y/sqrt(x**2 + y**2) - c*(x_dot + y_dot)*y/sqrt((H - x)**2 + y**2) + g*m - 1.0*k*(-L_1_init + sqrt(x**2 + y**2))*y/sqrt(x**2 + y**2) - 1.0*k*(-L_2_init + sqrt((H - x)**2 + y**2))*y/sqrt((H - x)**2 + y**2))*sin(beta)/2 - c_beta*beta_dot - 1.0*k_beta*beta - m*(-D*(-sin(beta)*beta_dot*x_dot + cos(beta)*beta_dot*y_dot)/2 + D*sin(beta)*beta_dot*x_dot/2 - D*cos(beta)*beta_dot*y_dot/2)/2 + m*(-D*(-sin(beta)*x_dot + cos(beta)*y_dot)*beta_dot/2 + D*sin(beta)*beta_dot*x_dot/2 - D*cos(beta)*beta_dot*y_dot/2)/2)*cos(beta)/(2*(-D**2*m*sin(beta)**2/4 - D**2*m*cos(beta)**2/4 + D**2*m/4 + Izz)) + c*(H - x)*(x_dot + y_dot)/sqrt((H - x)**2 + y**2) - c*(x_dot + y_dot)*x/sqrt(x**2 + y**2) - 1.0*k*(-H + x)*(-L_2_init + sqrt((H - x)**2 + y**2))/sqrt((H - x)**2 + y**2) - 1.0*k*(-L_1_init + sqrt(x**2 + y**2))*x/sqrt(x**2 + y**2))/m,\n",
    "              (D*m*cos(beta)*beta_dot**2/2 + D*m*(-D*g*m*sin(beta)/2 + D*(-D*m*sin(beta)*beta_dot**2/2 + c*(H - x)*(x_dot + y_dot)/sqrt((H - x)**2 + y**2) - c*(x_dot + y_dot)*x/sqrt(x**2 + y**2) - 1.0*k*(-H + x)*(-L_2_init + sqrt((H - x)**2 + y**2))/sqrt((H - x)**2 + y**2) - 1.0*k*(-L_1_init + sqrt(x**2 + y**2))*x/sqrt(x**2 + y**2))*cos(beta)/2 + D*(D*m*cos(beta)*beta_dot**2/2 - c*(x_dot + y_dot)*y/sqrt(x**2 + y**2) - c*(x_dot + y_dot)*y/sqrt((H - x)**2 + y**2) + g*m - 1.0*k*(-L_1_init + sqrt(x**2 + y**2))*y/sqrt(x**2 + y**2) - 1.0*k*(-L_2_init + sqrt((H - x)**2 + y**2))*y/sqrt((H - x)**2 + y**2))*sin(beta)/2 - c_beta*beta_dot - 1.0*k_beta*beta - m*(-D*(-sin(beta)*beta_dot*x_dot + cos(beta)*beta_dot*y_dot)/2 + D*sin(beta)*beta_dot*x_dot/2 - D*cos(beta)*beta_dot*y_dot/2)/2 + m*(-D*(-sin(beta)*x_dot + cos(beta)*y_dot)*beta_dot/2 + D*sin(beta)*beta_dot*x_dot/2 - D*cos(beta)*beta_dot*y_dot/2)/2)*sin(beta)/(2*(-D**2*m*sin(beta)**2/4 - D**2*m*cos(beta)**2/4 + D**2*m/4 + Izz)) - c*(x_dot + y_dot)*y/sqrt(x**2 + y**2) - c*(x_dot + y_dot)*y/sqrt((H - x)**2 + y**2) + g*m - 1.0*k*(-L_1_init + sqrt(x**2 + y**2))*y/sqrt(x**2 + y**2) - 1.0*k*(-L_2_init + sqrt((H - x)**2 + y**2))*y/sqrt((H - x)**2 + y**2))/m,\n",
    "              (-D*g*m*sin(beta)/2 + D*(-D*m*sin(beta)*beta_dot**2/2 + c*(H - x)*(x_dot + y_dot)/sqrt((H - x)**2 + y**2) - c*(x_dot + y_dot)*x/sqrt(x**2 + y**2) - 1.0*k*(-H + x)*(-L_2_init + sqrt((H - x)**2 + y**2))/sqrt((H - x)**2 + y**2) - 1.0*k*(-L_1_init + sqrt(x**2 + y**2))*x/sqrt(x**2 + y**2))*cos(beta)/2 + D*(D*m*cos(beta)*beta_dot**2/2 - c*(x_dot + y_dot)*y/sqrt(x**2 + y**2) - c*(x_dot + y_dot)*y/sqrt((H - x)**2 + y**2) + g*m - 1.0*k*(-L_1_init + sqrt(x**2 + y**2))*y/sqrt(x**2 + y**2) - 1.0*k*(-L_2_init + sqrt((H - x)**2 + y**2))*y/sqrt((H - x)**2 + y**2))*sin(beta)/2 - c_beta*beta_dot - 1.0*k_beta*beta - m*(-D*(-sin(beta)*beta_dot*x_dot + cos(beta)*beta_dot*y_dot)/2 + D*sin(beta)*beta_dot*x_dot/2 - D*cos(beta)*beta_dot*y_dot/2)/2 + m*(-D*(-sin(beta)*x_dot + cos(beta)*y_dot)*beta_dot/2 + D*sin(beta)*beta_dot*x_dot/2 - D*cos(beta)*beta_dot*y_dot/2)/2)/(-D**2*m*sin(beta)**2/4 - D**2*m*cos(beta)**2/4 + D**2*m/4 + Izz)\n",
    "             ]\n",
    "    return sysODE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Set up simulation parameters \n",
    "m = 1.0\n",
    "k = 10.0 \n",
    "# L = 1.0 \n",
    "g = 9.81 \n",
    "H = 20.0 \n",
    "c = 1.0 \n",
    "D = 1.0 \n",
    "Izz = 0.3333\n",
    "k_beta = 1.0\n",
    "c_beta = 1.0\n",
    "L_1_init = 18.6249\n",
    "L_2_init = 19.7227 \n",
    "\n",
    "\n",
    "# ODE solver parameters\n",
    "abserr = 1.0e-9\n",
    "relerr = 1.0e-9\n",
    "max_step = 0.01\n",
    "stoptime = 20.0\n",
    "numpoints = 100000\n",
    "\n",
    "# Create the time samples for the output of the ODE solver\n",
    "t = np.linspace(0.0, stoptime, numpoints)\n",
    "\n",
    "# Initial conditions\n",
    "x_init = 9.0                        # initial position\n",
    "x_dot_init = 0.0                    # initial velocity\n",
    "y_init = 17.0                   # initial angle\n",
    "y_dot_init = 0.0                # initial angular velocity\n",
    "beta_init = 0.0                        # initial position\n",
    "beta_dot_init = 0.0 \n",
    "\n",
    "# wf = np.np.sqrt(k / m1)                # forcing function frequency\n",
    "\n",
    "# Pack the parameters and initial conditions into arrays \n",
    "p = [m, k, g, H, c, D, t, Izz, k_beta, c_beta, L_1_init, L_2_init]\n",
    "x0 = [x_init, y_init, beta_init, x_dot_init, y_dot_init, beta_dot_init]\n",
    "resp = odeint(eq_of_motion, x0, t, args=(p,), atol=abserr, rtol=relerr,  hmax=max_step)\n",
    "freq, mag = CRAWLAB_fft(resp[:,2],t,False)\n",
    "np.savetxt(\"response.csv\", resp, delimiter=\",\")\n",
    "np.savetxt(\"sim_time.csv\", t, delimiter=\",\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.35000350007\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ 17.        ,  17.        ,  17.        , ...,  17.00003265,\n",
       "        17.00003265,  17.00003265])"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEpCAYAAABSuB9PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHoJJREFUeJzt3XmUZEW1qPGvGpBBGkEspqdC47BldMQBUQSVe1GUi8Nz\nAK6A0qgoIIOAICiKAjIIyCwCasvlIoqgIjgh+EREQQSFrVelUdurzdwIMnW9P+KUnZ3UlKdPVWZ2\nfb+1clVlnGlXrKrcFSfiRAwMDQ0hSVIdM7odgCSpf5lEJEm1mUQkSbWZRCRJtZlEJEm1mUQkSbUt\n2+0ApImKiHOAd7UVPwz8FbgUOCwz76l57pWBFTLzjiWLEiLiXcA51ds3Z+bXR9nvs8BewG2ZuX6N\n6ywHDGbmvLbrviozr6oVvNQhWyLqN0PA3sBO1WtP4GvAu4FvR8RApyeMiBcAtwIbNhgnlFjfOMb2\n7at9OhYRTwduAl7TUnwVpU5uqXNOqQ5bIupH38jM21sLIuJ3wCnAtsC3OzzfJsDaDcXW6o/A6yNi\nIDMXSxYR8XxgXeDvNc89C3h2a0Fm/rG6pjRlbIloafFDYADYqMaxHbdeJmAI+AawOrDFCNt3AH5H\n/VbDZMQsdcyWiJYWT6++/r61MCI2AD4FvAp4AnADcERmXlFtPxw4nPKhf2VE/Kt/IiK2BvYHXgys\nQmk1fBM4MDPvnUBMVwB7UG5pXd22bQfg68BL2g+qblUdCfwbMBNI4HOZ+flq+3DfxxBwbkSck5nL\nRMQuwBdo6ROJiBWBw4C3A+sA84D/qurgwWqf4eOeBxwE/DuwHPA9YJ/MnDuBn1XTlC0R9aMnR8Tq\n1WvtiHgNcCLwc+CS4Z0iYhPgGuA5lA/lj1D+cfp2RLy12u0i4Mzq+yOBfapjt6EkgZWAjwIfBK4F\nZgNnTDDOB6tzbN9aGBHPovS/PK7DPSLWq36ON1TX2R+4EzgzIo6qdruKkhgHqn12qsqHaOljqTre\nvwccAHyX0on/Q+BA4PKIWKblOCh19yTgYOA0YDvgggn+rJqmbImo3wwA149Q/gCwVWY+2lJ2MqX1\n8PzM/CdARJxM+SA9MSK+npk3R8Q1wO7Ad1tGNe0DzAVenZmPVWVnRMRPKP+pT9TFwPYRsUFmDt+6\n2gGYl5nXRkT7/kcBqwEvyswbq7JTIuISYP+IOC8zb4mI71KS4jWZef4o13438DJg78w8ueVn+A1w\nTPUzn96y/88y8/8Ov6lGrO0REc/IzMVaeNIwWyLqN0PAOymjkl4DvA54H6VD+erqFhQR8WTglZRO\n9icOt1woH9AXA2sCm41xnddTPsiHEwjV8fcBK3cQ7zeBhSzeGvmPKobFRMSM6ue5vCWBDDuS8vc6\n1mivdm8A7gVObSs/kfJztMY0BFzYtt8vq69rdXBNTTO2RNSPfjLC6KwLgf+htD42Ap5Rbfog5TZO\nuyFKP8o1I10gM4ci4plVf8Hw+f5Py7ETkpl3RcSPKR/+R0XEOpR+kENG2P0plASVI2wbbsWsO9Fr\nU0Zw/aE1EVYxPRIRfxjhXPPb3j9UfV0GaRS2RLRUyMy7KLepnhMRT2LRB98pLGq1tL5eC1w52vki\nYn9KH8hWlA/1oykf/nNqhHcx8OKIWJPy3//dwI9G2G+sEVfDf6sPd3Dd8c7Xfq6FHZxbAmyJaOky\n/EG7ELit+v7RzPxB607ViK1ZlH6Ux4mI5YGPAd8Html9xqNKBJ26GDiecntpe+DSzBzpA3s+8A/K\nQIB2w2V/6uC6twEvjYhl2m7LLUf5+X2qXUvMloiWCtWH+9bADZm5IDP/lzLKaZeIWLtlv2Upw2Mv\nZNE/UcMfsMN/DytSRmX9ri2BPI/SzzLcfzEh1RDZXwHvALakPGE/0n4LgcuAbaprtTqQkhy/NUrM\nI7mUMtpqz7byPSm3zS6d4I8gjcqWiPrRDhExPMfVAKVvY3fKB/9HWvbbi9Ka+EVEnEoZKvtOSof6\nQZl5d7Xf/Oo874+ItTPz/Ii4FtgtIhZQbmdtQhnt9BjlGYqZlE7r0bTfSrqY8jzK/ZRhv6M5iHIL\n7UfVSLK/Am+iPOdyXGbe2hIzwM5VQjtvhOt+njLX2PERsSklqW4G7ELpCzp7jHjHK5cAWyLqT8cD\nX6xeX6Aki5uA1ww/RAiQmT8FXg5cB+xLGda6IvCuzPxMy/m+T3ke4nXAyRHxBOCtlCfOdwVOAF5N\neTZjx+qYrceJsb3z/eKq7LLMfGi0fTPzD5S+l29RHlQ8mvKg426Z+eGW/RI4CXhhFd/TRzjXw1Wc\nx1P6gU6gtKQ+yeJDl0eKd7xyCYCBoSF/RyRJ9dgSkSTVZhKRJNVmEpEk1WYSkSTVtlQN8X300ceG\n7r57xOfHVMNqq62E9dkM67JZ1mezBgdn1h7KvVS1RJZd1il+mmR9Nse6bJb12TuWqiQiSZpaJhFJ\nUm0mEUlSbSYRSVJtJhFJUm0mEUlSbSYRSVJtJhFJUm09+cR6RKxEWUfhTZSFhq4B9svMW7oamCRp\nMb3aEjmJspjOm4GXAv8ELqsWC5Ik9YheTSLbA6dm5k+rFdwOoazctmF3w5IkterJ21mU9aPfFhH/\nTVnH+j2U9bH/0NWoJEmL6dUkMhv4MvA34DHgH8A2mXlfV6OSJC2mV29nPQv4K7AtsDlwOXBRRKzT\n1agkSYvpuSQSEesBZwJ7ZeblmXkdsCOlc/1D3YxNkrS4Xryd9SJKcvvFcEFmPhoRNwDPHO/gwcGZ\nkxja9GN9Nse6bJb12Rt6MYn8ufq6KfDLlvINgW+Pd/D8+QsmI6ZpaXBwpvXZEOuyWdZns5YkIfdi\nEvkZcC1wbkTsCdxBuY31NODkbgYmSVpcz/WJZOZCYDtKIjmf8rT6+sAWmfmnbsYmSVpcL7ZEyMy7\ngD26HYckaWw91xKRJPUPk4gkqTaTiCSpNpOIJKk2k4gkqTaTiCSpNpOIJKk2k4gkqTaTiCSpNpOI\nJKm2npz2pK7LfzqX++//Z7fDWGqsvPIK1mdDrMtmWZ/NGRiAN706ah+/VCWRz134y/F3kiQtZkmS\niLezJEm1DQwNDXU7hsbccc+DQ3feeX+3w1hqrL76ylifzbAum2V9NmfGjAGeNespA3WPX6puZz1l\n1RUZeuTRboex1LA+m2NdNsv67B3ezpIk1WYSkSTVZhKRJNVmEpEk1WYSkSTVZhKRJNVmEpEk1WYS\nkSTV1vWHDSPidGBGZs5uKdsGOBoI4LfAQZn5nS6FKEkaRVdbIhFxBDC7rWxD4BvABcDzgEuAiyNi\ng6mPUJI0lq60RCJiFnA2sBEwt23zXsA1mXlU9f6wiNgC2Bt479RFKUkaT7daIpsDtwObALe1bXsF\ncGVb2ZVVuSSph3SlJZKZc4A5ABGPm8f+qcBf2srmAU+b/MgkSZ3oxdFZKwHtS5Y9BKzQhVgkSWPo\n+uisETwILN9Wtjzwj4kcPDg4s/GApjPrsznWZbOsz97Qi0nkT8DabWXr8PhbXCOaP39B4wFNV4OD\nM63PhliXzbI+m7UkCbkXb2f9GNiyrWwr4KouxCJJGkMvtkROBn4eER8Dzgd2BF6Mw3slqef0Qktk\nsUXeM/NmYAfgzcANwHbAdpmZXYhNkjSGgaGhofH36h9D3idtjvedm2NdNsv6bNbg4MyBusf2QktE\nktSnTCKSpNpMIpKk2kwikqTaTCKSpNpMIpKk2kwikqTaTCKSpNpMIpKk2kwikqTaTCKSpNpMIpKk\n2kwikqTaTCKSpNpMIpKk2kwikqTaTCKSpNpMIpKk2kwikqTaTCKSpNpMIpKk2kwikqTalu12ABFx\nOjAjM2e3lH0A2BN4GnAbcEJmnt2dCCVJo+lqSyQijgBmt5W9D/g0cASwCXACcGpE7Dj1EUqSxtKV\nlkhEzALOBjYC5rZt3gM4OTPPr96fHREvA3YF5kxdlJKk8XSrJbI5cDulpXFb27YPAme0lS0EVpv8\nsCRJnehKSyQz51C1KiKifdvVre8j4unAO4ATpyo+SdLEdJREImJn4I7MvCwiNgW+ROn8vhD4YGY+\n3GRwETEIfAuYBxzd5LklSUtuwkkkIvYDjgEOAy4DTgPWAM4FdgHuBg5qKrCIWL+6zvLAlpm5YCLH\nDQ7ObCoEYX02ybpslvXZGzppibwHOCYzj4yI9YCXAXtm5mkRcStwCA0lkYh4ASWBzAe2ysx5Ez12\n/vwJ5RpNwODgTOuzIdZls6zPZi1JQu6kY30W5YMd4PXAEHBp9f5WSqtkiUXEc4ArgN8DW3SSQCRJ\nU6uTlsh8FiWK1wG3Zuafq/ebAv/bUExfBB4E/hNYPiLWrMofzcw7G7qGJKkBnSSRS4GjI+I1wLaU\n21dExL6UfpK6T5QPDX8TEc8CXli9zbb9fg88u+Y1JEmToJMksi9lmO2WwOnAsVX5bOAS4NA6AWTm\n1i3f/w5Yps55JElTb8JJJDP/SXmavN1zM/Oh5kKSJPWLjh82rKYgeS2wDvApYIOIuCEz/950cJKk\n3tbJcyJPoDxl/mbgYWA54EzgAGDDiHhFZv5+UqKUJPWkTob4fhLYBtgeWBUYqMrfA9wLHNlsaJKk\nXtdJEtkRODgzLwUeGS7MzNuAjwOvajQySVLP6ySJPBn4n1G23QGssuThSJL6SSdJ5NfA20fZti3w\nmyUPR5LUTzoZnXUkcFFEPJny4OEQ8PKI2An4ALDTJMQnSephE26JZObXKYniBcBZlI71Eyl9JXtm\n5n9PSoSSpJ7V0XMimfkV4CtRVpJanTIq65bMXDgZwUmSelvHy+NGxArAmsDTKYtFrdN0UJKk/tBR\nEomIPSmJ40rKg4ezgDMi4nsR8cTmw5Mk9bIJJ5GI2A04ibKS4atZ9LDh2cBmlGdFJEnTSCctkQOA\n4zJzX+Cq4cLM/BplBt+3NBybJKnHdbqy4RWjbLsJWGvJw5Ek9ZNOksifgRePsu351XZJ0jTSyRDf\nLwAfjYgHgG9WZStGxBspqxye1HRwkqTe1kkS+TSwLnBc9YJFfSP/hbP4StK008nKhkPAHhFxLLA1\nix42vCozb5qk+CRJPayTRamuBw7NzG8Dv5u8kCRJ/aKTjvVnAg9MViCSpP7TSRI5H/hQRKwxWcFI\nkvpLJx3r6wFbAX+NiL8B97dtH8rM6DSAiDgdmJGZs0fYtixwHXBDZu7W6bklSZOrkyQyjzJfVmMi\n4ghgNvD5UXb5BPBc4IYmrytJakYno7N2beqiETGLMufWRsDcUfZ5ObAr8KumritJalYno7OePsbm\nhcD9mXnPBE+3OXA7ZbndC0a41hOB8ygrJu4x0RglSVOrk9tZt1GWxB1VRNwFnJSZnxhrv8ycQ3Vr\nrKxv9TgnAtdm5lcjwiQiST2qk9FZuwAPA9+h3GbaFngXcDGlJfIJ4Bzg4GrdkVqqaVT+HXh/3XNI\nkqZGJy2RdwJfzszd28q/HBEnAy/KzO0i4k7gfcApnQYTEYPAmcCumXlvp8dLkqZWJ0nklcD2o2z7\nBnBJ9f01wOE149kWGAQuiIjhRa9WAIYi4i2Zucp4JxgcnFnz0hqJ9dkc67JZ1mdv6CSJ/J3SIf7d\nEbZtDtxZfb8qcF/NeC4CftxW9kXgr8CHJ3KC+fMX1Ly02g0OzrQ+G2JdNsv6bNaSJOROksgZwMcj\nYkXga8B8YA1K62Q/4MiIWAs4mJaVDzuRmf8A/tBaFhEPAgsy8491zilJmjydPCfy6YhYAdifslQu\nlHXW7wM+Wb12orREduwghjFHfEmSetfA0FBnn+ERsRLwMkrfxV+AX2bmgmrbspn5aONRTtyQTdzm\neMugOdZls6zPZg0OzhwYf6+RdTLEd9hC4JHq+5uBJw1v6HICkSRNsY6SSPX8xzzgSsrDgrOAMyLi\ne9VT5pKkaWTCSSQidqOso34u8GpKfwiUObA2Az7edHCSpN7WSUvkAOC4zNyXltFXmfk14FDgLQ3H\nJknqcZ0kkVnAFaNsuwlYa8nDkST1k06SyJ+BF4+y7fnVdknSNNLJw4ZfAD4aEQ8A36zKVqwmTDyE\n0l8iSZpGOkkinwbWBY6rXlD6RgYo668f2WxokqRe18kT60PAHhFxLLA1sDpwL3BVZt40SfFJknrY\nhJJItVTtbGALFnWg/wW4Grh2ckKTJPW6cTvWI+IEym2rt1CWtL20ev0NeAdwbUR4K0uSpqExWyIR\nsSuwN3A0cGRm3t+2fRXKrL0HRcQNmfnVSYtUktRzxrudtTtwXmYePNLGzLyPshzu2pTVDE0ikjSN\njHc7a0PK2iHj+TrwvCUPR5LUT8ZLIiuzaMXCscynZTZfSdL0MF4SmQFMZHr3x1g0IaMkaZqYyLQn\nrjwoSRrRRJ4TOTki7htnn1WaCEaS1F/GSyJXUVoiy42z34O0TA8vSZoexkwimfmqKYpDktSH6qyx\nLkkSYBKRJC0Bk4gkqbZO1hOZFBFxOjAjM2e3lG0AfJYya/BdwDmZeViXQpQkjaKrLZGIOIIyxXxr\n2erAlcAdlKlU3g98MCL2nfIAJUlj6kpLJCJmAWcDGwFz2zbvRVnsaufMXAj8LiKOAzYHjp/SQCVJ\nY+rW7azNKWuTvB24oG3bNsDXqwQCQGZ+cgpjkyRNUFeSSGbOAeYARET75mcDF0bEScCbgAXAecAx\nrYlFktR9vTg6axXgEOARYDvgk8CBgB3rktRjuj46awSPADdm5n7V+19GxFrAocDHxjt4cHDmJIY2\n/VifzbEum2V99oZeTCJ/AW5uK/sNsEpErJaZd4918Pz5CyYtsOlmcHCm9dkQ67JZ1mezliQh9+Lt\nrKuBzdrKNgHuGi+BSJKmVi+2RI4FrquG9Z4KPBc4CDiuq1FJkh6nF1oiiy16lZm/AV4DvAS4CTiB\nMjLryC7EJkkaQ9dbIpm59Qhl11CmPJEk9bBeaIlIkvqUSUSSVJtJRJJUm0lEklSbSUSSVJtJRJJU\nm0lEklSbSUSSVJtJRJJUm0lEklSbSUSSVJtJRJJUm0lEklSbSUSSVJtJRJJUm0lEklSbSUSSVJtJ\nRJJUm0lEklSbSUSSVJtJRJJUm0lEklTbst0OICJOB2Zk5uyWsrcCHwXWB24Djs3Mc7sSoCRpVF1t\niUTEEcDstrJXAHOAk4CNq69nRcS2Ux+hJGksXWmJRMQs4GxgI2Bu2+Y3Ajdm5uer92dGxLuBfwMu\nm7ooJUnj6VZLZHPgdmATyu2qVvOBjSLiVQAR8UpKsrluCuOTJE1AV1oimTmHcsuKiGjffAolyfwg\nIh6jJLpjq2MkST2kF0dnrQGsBewPvBDYC/hAROzSzaAkSY/X9dFZIzgLuD4zj6/e/yoi1gCOAc4d\n7+DBwZmTGNr0Y302x7pslvXZG3oxibwU+Epb2bXAoRHxpMy8d6yD589fMGmBTTeDgzOtz4ZYl82y\nPpu1JAm5F29n/RnYtK1sE+DO8RKIJGlq9WJL5ETg+Ii4Bbic0sl+MPCxbgYlSXq8XkgiQ61vMvO0\niHgI2Bs4lvIcyUGZeXo3gpMkjW5gaGho/L36x5D3SZvjfefmWJfNsj6bNTg4c6Dusb3YJyJJ6hMm\nEUlSbSYRSVJtJhFJUm0mEUlSbSYRSVJtJhFJUm0mEUlSbSYRSVJtJhFJUm0mEUlSbSYRSVJtJhFJ\nUm0mEUlSbSYRSVJtJhFJUm0mEUlSbSYRSVJtJhFJUm0mEUlSbSYRSVJtJhFJUm3LduOiEbEG8Bng\ntcCKwLXAfpn562r7NsDRQAC/BQ7KzO90I1ZJ0uimvCUSEQPAxcAzgTcALwPuBb4fEatFxIbAN4AL\ngOcBlwAXR8QGUx2rJGls3WiJPBd4CbBBZv4WICJ2Bu4CXg9sAVyTmUdV+x8WEVsAewPv7UK8kqRR\ndKNP5HZgu+EEUllYfV0NeAVwZdsxV1blkqQeMuUtkcy8C7isrXhvYAXgCuCTwF/ats8Dnjb50UmS\nOtH10VkR8UbgU8BxmZnASsA/23Z7iJJkJEk9pKtJJCJ2Ab4KnJ+ZB1bFDwLLt+26PPCPKQxNkjQB\nXRniCxARhwCfAE7KzH1aNv0JWLtt93V4/C2ukQwMDs5sKEIBWJ/NsS6bZX32hq60RCLiw8ARwKFt\nCQTgx8CWbWVbAVdNRWySpIkbGBoamtILRsSmwC+Ac4FD2zYvANYHfg4cBZwP7AjsB7yg6jORJPWI\nbrRE3lZddzfKqKvW1z6ZeTOwA/Bm4AZgO8qQYBOIJPWYKW+JSJKWHl0f4itJ6l8mEUlSbV0b4tuU\niJgBHAm8C5gJfAfYMzP/3tXA+lQ10eWvgSFgoCoeAl6RmT/pWmB9KCJOB2Zk5uyWMmeormGUuvwZ\n8KKW3YaAs1v30SKTNXv60tAS+TiwM7ATZX6tp1IeYFQ9mwDzgbVaXmtTfuE0QRFxBDC7rcwZqmsY\nqS4rGwLvYPHf032nMLS+MZmzp/d1SyQilgP2Aj6QmT+oyt4O/DEiXpqZP+1qgP1pY+A3mTm/24H0\no4iYBZwNbATMbdu8F85QPWFj1WVErE/5b/qn3nWYkEmbPb3fWyLPA1YGfjRckJlzgdtw1t+6NgZu\n6XYQfWxzykzVm1B+D1s5Q3VnxqrLjYEHq793jW/SZk/v65YI5dYVOOtvkzYGVoiIa4D1gJuBj2Tm\ndV2Nqk9k5hxgDkBEtG9+Kv6uTtg4dbkxcG9EfIUyw8WdwDnAZzPT5xbaTObs6f3eElkJWJiZj7WV\nO+tvDRGxAmXGgJnA/pR7p/OAH8UIf8XqmDNUN2cj4ImUD8ZtgM9R+kcP62ZQ/aLJ2dP7vSXyIDAj\nImZk5sKWcmf9rSEz/xkRqwIPZeYj8K+Zll8IvJ/yn4vqc4bq5uwMrJyZ91Xvf1397n6Ekkw0iupv\n+kzgK03Mnt7vLZE/VV/rzvqrNpl5/3ACqd4PUYb8estlyS3JDNVqkZkLWxLIsJuAmRGxSjdi6gfV\n7OlfAE7NzF1aNtX+3ez3JHIjcD8ts/5GxHqUe/nO+tuhiHhBRNwbEc9vKZtBGcBwc/ciW2o4Q3VD\nIuKaiPhsW/FmwLwRkouYvNnT+/p2VmY+HBGnAsdGxJ2U5xtOAX6YmT/rbnR96Ubgj8AZEfEBSlP2\nQGB14KRuBraUOBn4eUR8jEUzVL8Yh/fW8TXg4xHxC+D/UT7wDqAMo1abavb0IymtkLMjYs2WzQtY\ngt/Nfm+JQJlOfg7wJeD7lA/Bt3Y1oj5VDVDYFkjKw0Y/BdagPK1+Rzdj61OLjRJyhuol0l6Xn6H0\nfxxCaSUfQJkF/JwuxNYPJm32dGfxlSTVtjS0RCRJXWISkSTVZhKRJNVmEpEk1WYSkSTVZhKRJNVm\nEpEk1WYS0bQWEa+PiPNa3m8ZEQsjYvMux/WRiDhrgvuuGhF/jIh1JzsuqZ1JRNPdPiw+ueQvgJdS\npoDpiojYGPgQcNBE9s/MeyhrZ587iWFJIzKJSC2qWYx/lpndnJ79KOC8zLyzg2POAjaMiP+YpJik\nETntiaatiPghi2YuHaJM4jcA/BDYIjN/EhGHA2+hrFHxCcoM0TcAuwDPBj5NWcjrV8B7M/PGlvPv\nQJnbaSPKWtZzKDOoPjxGTBtX53pRZl5flQ1U134nZXrueZRJ8g7PzEdbjj0VeGFmvqR2pUgdsiWi\n6ex9wHXA9ZRbWNdX5e3/Wa1HWT70EMrsps8CvgUcR/lwfxuwLvDF4QMi4p3ARZSEsD1lBtU9qJZ7\nHcOOwNzhBFI5iDKb6uHAa4FTgQ8DB7cd+1XgRRHxjHGuITWmr6eCl5ZEZt4aEfcBywyvIT/KKsAr\nAbtn5tXVPq8C9gS2zswfVWXHAp+JiJUy8wHKLalLMnO36hxXRMSfgIsj4mWZec0oYW1FSWytXgn8\nPDO/VL2/OiIeAO5p2+/nlJbUVsDvx68BacmZRKSJaV2f5m8jlA33X6waEU8DnkpZ72KZln2+CzxC\naU2MlkTWB37UVvZD4KiIuIoyRf+3MvPU9gMz876IuIfScpKmhLezpPE9lpkPtRdm5oOj7L969fVM\nStIYfj0ALEfp1xjNk2hb1zozj6G0fFaktHB+HRE3VS2idv+oziFNCVsiUvPurb7uQ1l1r91YC3zd\nAazaXpiZpwGnRcRTgNdRFmP7akSsWS0mNmy1cc4vNcqWiKa7x8bfpWO3UJZqnpWZ1w+/KCO0jgE2\nGOPYuZRbYf8SEVcNryeemXdk5heBz1ESxhNb9luV0n9ze5M/jDQWWyKa7u4BtoiIrShDd6F0TteW\nmQsj4lDglIgYAi4DnkIZXfWkluuM5Apg97ayK4EDI+JvwE8oSWY/4HuZeV/LfltQRpZdsSTxS52w\nJaLp7lRKf8W3KR3e8PghviM9TDXmA1aZeRawE+U5lEuAz1Kegn9FZv59jEMvAtaKiOe2lB1OGWK8\nKyUhHVvF+7a2Y7cFfpaZfxkrNqlJPmwo9ZiIuBS4PTP37OCYFSkPIf5nZl46acFJbWyJSL3nUODt\nEbFmB8e8F/i1CURTzSQi9Zhq6pTjgU9NZP+IWA3YG9h5MuOSRuLtLElSbbZEJEm1mUQkSbWZRCRJ\ntZlEJEm1mUQkSbWZRCRJtf1/WkyZKeimQfYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f7471d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcMAAAEpCAYAAADvWUk6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcXFWd9/FPdXfS2RMCnbCERVR+KKKgOJKQsGSIIDLj\n4OiIIgITgQdZFAHFBQVZXEAFQUScqA9oUFFRHB4R2YJhEQKyufxQgRAIkIQmWyfp9FLPH+fcrtvV\nVV1VSXeq0/f7fr361VX3nDrn1CXkl3PuWXL5fB4REZEsa6h3A0REROpNwVBERDJPwVBERDJPwVBE\nRDJPwVBERDJPwVBERDKvqd4NEJH6M7MvAl/sJ0se2Bc4qsp8ZwLHVVH1D939v6ttp8hgUTAUkUQe\nuBj4W5n0xTXkuwb4feraLOAk4FrgD6nr/9zUxooMJAVDEUm73d3vKZdoZlXlA/4Yf5LPjSAEw/vd\nff5ANFRkIOmZoYiIZJ6CoYiIZJ6GSUUkbaKZbVvi+ip379yEfCJbBQVDEUnkgF+XuJ4HDgHuqTGf\nyFZDwVBEEnngLODxEmmPbUI+ka2GgqGIpD1SYZZorflEtgqaQCMiIpmnYCgiIpmnYCgiIplXl2eG\nZtZA2M7pOGA8cCtwqrsvK5N/P+Bywp6HzwMXufv1qfTRwBWEfRObgBuBM929LZXnGOA8YBfCQ/7T\n3X1RKv21wFXATKAVuNLdL6ulzWZ2LmGXje2BPwNfdPf/l0pvAb4NzAE2Aj8APuvu3dXeO5GtWK7e\nDRApp149wwuAY4EPE/YsnAb8vFRGM9uOEHgWEYLhlcA8Mzs0le1aYAZwBHAkcDBhb8SkjEOBecCl\nsYwngNuSdVJxq6hbgVXA24FPA+eb2dwKbf5Fqo6zgE8BZwBvBH4J/NrM9kmV8UtgSvz8ccAJsVyR\nLMjXuwEi5eTy+S375zMGnhXAaUnvzsx2BZ4BZrj7A0X5PwPMdffXpa59H9jR3Q83s2nAs8Ah7v6H\nmH4gcBcwzd1fNLNbgaXJ7vhmlgP+DvyPu3/FzD4IfBeY6u7rY54vAB9y9z2rabOZfQ54wd1/mGpn\nK3CBu19hZtOBhcBr3P25mP4R4FtAi7t3DMT9FRGR2tWjZ7gPMA5YkFxw98WEgDarRP6Z9F3Eezdw\nQHw9A+gC7kul3xuvzYyB74D4maS+fCwzqW8msCgJhKk6Xh+HNvet1GZ3vzgJhGbWbGYfA0YTgnJS\nx+IkEKbqmBDviYiI1Ek9guG0+PuFoutLgZ3L5C+Vd4yZTQZ2Apa5e1eSGF8vi+VNAsZWqK9cHcQ8\nO1XbZjM7ClhHGM69yN2ThcmV6hARkTqpRzAcA3Sng1fUDowqk39DibzE/KXS0+WNie9LlZHUV66O\nXKqMatv8AKGn90ngi2b20XJ1xD0c8yXKEBGRLagewXA90BBnZ6Y1A21l8jeXyEvMXyo9Xd76os+U\nqq9cHflUGVW12d1fdPcn3P0K4H+As8vVYWZNhIBb6nuLiMgWUo9guCT+3qHo+o70HUZM8pfKu9bd\nV8X0KfHZIABm1kiYtfm8u7cSgk1/9ZWrg5inYpvN7Ii4PCPtCQpDrJXqKCsfZjnpRz/60Y9+avup\nWj3WGT4GrAUOAuYDmNluwG6U3u1+IXB80bXZhEkyxN9NwHQKk2hmEXpcSZ77Yn0/jvXlgAMJM0iT\nOj5kZqPcPRnKnA24u68ws9X9tDmZVHMpcAdhaUXiHcBfUnV8xcx2cvck+M0GVgOPlvjePXK5HMuX\nr+kvS2a0tIzXvYh0Lwp0Lwp0LwpaWsZXnXeLB0N332hmVwOXmdkrwHLCQvS73P3BuIxhMtAalxvM\nA84xs+8QFtbPAY4GDovlLTWzGwlrD+cServXAte5+0ux2m8AN5vZo8CdhB33J8SyAW4CLgLmm9l5\nwJsJw5unVNHmh1J1fNvMHiYEvvcCHwT+PZZxv5k9APzUzE4nLMz/KvB1nf8mIlJf9Vp0/3lCL+16\nQm/qGeD9MW0GYZbldIC4w8vhhOUNjwAfA4519wWp8uYSen+3EALb7TEfsYzfEXaG+STwMLAnMCcO\noRJ7g4cTAuSDwCXAueldbiq0GXefR+gVfoYwPPoh4L2x7sRRwMuEHvA84Fp3v7DquyYiIoNiiy+6\nl02S17BHoCGgAt2LAt2LAt2LgpaW8VVvAaiNukVEJPMUDEVEJPMUDEVEJPMUDEVEJPMUDEVEJPMU\nDEVEJPMUDEVEJPMUDEVEJPMUDEVEJPMUDEVEJPMUDEVEJPMUDEVEJPMUDEVEJPMUDEVEJPMUDEVE\nJPMUDEVEJPMUDEVEJPMUDDMun8/XuwkiInWnYJhhK1au55zv3MeVv3hcQVFEMq2p3g2Q+lnky2ld\n3U7r6nZWr+tg4tiR9W6SiEhdqGeYYSvXtve87ujsqmNLRETqS8Eww1a1bex53dWlYVIRyS4Fwwxb\nuabQM+zs6q5jS0RE6kvBMMPSw6Sd6hmKSIYpGGbYytQwaWe3eoYikl11mU1qZg3AxcBxwHjgVuBU\nd19WJv9+wOXAvsDzwEXufn0qfTRwBXAU4TvdCJzp7m2pPMcA5wG7AI8Bp7v7olT6a4GrgJlAK3Cl\nu19WS5vN7DTgVGBn4Fngm+4+L5V+CvBtIA/k4uVOd9/i0zjXt3fSvrEwaUbPDEUky+rVM7wAOBb4\nMDALmAb8vFRGM9uOEHgWEYLhlcA8Mzs0le1aYAZwBHAkcDBwTaqMQ4F5wKWxjCeA28xs25g+Itax\nCng78GngfDObW6HNv0jVcQrwZeBLwN7AN4GrYxBO7A38Gtg+9bNT/7dqcKQnzwB06JmhiGTYFu8Z\nxsBzBnCau98Zrx0NPGNm+7v7A0UfORFY6e6fiO+fMrO3AmcDt5vZNOCDwCHu/lAs76PAXWb2KXd/\nMeadn/TSzOxkYHYs+yvA+4CpwAnuvh74m5ntAZxDCLzVtPlkQm/yhtjOeWY2HTgB+HG89ibgDndf\nvvl3cvOkJ88AdCkYikiG1aNnuA8wDliQXHD3xYRhxVkl8s8E7im6djdwQHw9A+gC7kul3xuvzTSz\nXMx7d6q+fCwzqW8msCgGwnQdrzezFkJvslKbTwe+W9TObmCb1Pu9gL+W+I5bXHryDGgCjYhkWz2e\nGU6Lv18our6U8KytVP5HSuQdY2aTCcOMy9y95wGYu3eZ2bJY3iRgbJn69kvVUSqdWEYylFm2ze7+\nh3SCme1C6LFeEd/vSAiMR5jZBbFNC4Ck97pFrVzbe5hUSytEJMvq0TMcA3Sng1fUDowqk39DibzE\n/KXS0+WNie9LlZHUV66OXKqMqtsce5O3EILlV+PlvQgTZ9qBDwDHA3sAd5hZc4n2D6rinqEm0IhI\nltUjGK4HGuLszLRmoK1M/uJgkbxvK5OeLm990WdK1VeujnyqjKrabGa7AwsJM07f6e5rANz990CL\nu5/s7o/HZ4/vAYww8WeL6jtMqp6hiGRXPYZJl8TfO9B72HFH+g5DJvl3KLq2I7DW3VeZ2RJgipnl\n4rNAzKwRmAI87+6tZtZWpoykviWEXlpxOjFPsvSh3zbHiT2/BZYTJvQsTeXF3VuL3r9kZisoPTzc\nS0vL+EpZarJuY+9O7qgxIwe8jsGytbRzS9C9KNC9KNC9qF09guFjwFrgIGA+gJntBuxG34kyEHpZ\nxxddm02YJEP83QRMpzCJZhZhiDPJc1+s78exvhxwIIUJLwuBD5nZKHdPhktnA+7uK8xsdT9tXhDf\n7wncBjwFHOHuK9MNNrPTgXOBXZLhVjPbFWgBnix1o9KWL19TKUtNlreu6/V+5cr1A17HYGhpGb9V\ntHNL0L0o0L0o0L0oqOUfBVs8GLr7RjO7GrjMzF4h9KK+Ddzl7g/GZQyTgVZ37yCsDzzHzL5DmIwy\nBzgaOCyWt9TMbiQsZZhLGPq9FrjO3V+K1X4DuNnMHgXuBM4CJsSyAW4CLgLmm9l5wJsJyzFOqaLN\nD8UyriMMp34EaDazqfF6p7u/QniGeFFs55eB7QgbCdyTLNfYklYWrTPUDjQikmX1WnT/eUIv7Xrg\nDuAZ4P0xbQZh4sl0gLjDy+GE5Q2PAB8DjnX3Bany5hJ6f7cQAtvtMR+xjN8BJwGfBB4G9gTmJMOW\nsTd4OCFAPghcApyb3uWmvzab2euBtxGGTT22P/m5P9bxNCGQ7wz8EfgV8CjhueEWVbz7DGhphYhk\nW04nnG8V8gM57PHiK2187nt/7HXt3w/Yjf+YtfuA1TFYNARUoHtRoHtRoHtR0NIyPlc5V6CNujNo\nVdEaQ9B2bCKSbQqGGVS8rAK0zlBEsk3BMIPSu880NoRRBK0zFJEsUzDMoKRnOLKpgbGjRwCaQCMi\n2VaX8wylvpJgOHHcSLq7QxDUqRUikmXqGWZQMoFm0rhmmhrDH4HObvUMRSS7FAwzKOkZ9gqGneoZ\nikh2KRhmUDKBZuK4kTQ2agKNiIiCYcasb++kvSPsPrONhklFRAAFw8xJrzGcOG4kTXFphSbQiEiW\nKRhmTHqN4aRxzTQmPUMtrRCRDFMwzJhVqZ7hpHHNjGhKgqF6hiKSXQqGGdO7ZzhSO9CIiKBgmDnp\n3WdGNzcVJtBomFREMkzBMGPSawxzuRxNWlohIqJgmDWr4wn3E8aNBOiZQNOlpRUikmEKhhmzMe40\nM2pEI0BqmFQ9QxHJLgXDjEm2XUuCYGGYVD1DEckuBcOMSXaaSYKgeoYiIgqGmZMEvb49QwVDEcku\nBcOMSYJeskF3U0P4I5DP03O2oYhI1igYZkzybHBE7BkmQTGkqXcoItmkYJgxXT09w2SYtPBHQJNo\nRCSrFAwzpqPnmWHvCTQAnd3qGYpINikYZkxXVzKbtPcEGtBp9yKSXQqGGdKdz/fsNNNUaphUE2hE\nJKOa6lGpmTUAFwPHAeOBW4FT3X1Zmfz7AZcD+wLPAxe5+/Wp9NHAFcBRhO90I3Cmu7el8hwDnAfs\nAjwGnO7ui1LprwWuAmYCrcCV7n5ZLW02s9OAU4GdgWeBb7r7vFR6C/BtYA6wEfgB8Fl33yJdsvQB\nvkmPMD2BRgf8ikhW1atneAFwLPBhYBYwDfh5qYxmth0h8CwiBMMrgXlmdmgq27XADOAI4EjgYOCa\nVBmHAvOAS2MZTwC3mdm2MX1ErGMV8Hbg08D5Zja3Qpt/karjFODLwJeAvYFvAlfHIJz4JTAlfv44\n4IRY7haRniDTGJdUJEsritNFRLJkiwfDGHjOAD7j7ne6+6PA0cBMM9u/xEdOBFa6+yfc/Sl3vwr4\nEXB2LG8a8EHgFHd/yN3vBT4KfMjMdohlnA3Md/d57u7AyYTe34kx/X3AVOAEd/+bu/8E+BpwToU2\nH5Bq88mE3uQN7v5M7BFeTwh4mNl0QsD+iLs/6e63xvJPj+UPuvTSieRQ396zSdUzFJFsqkfPcB9g\nHLAgueDuiwnDirNK5J8J3FN07W7ggPh6BtAF3JdKvzdem2lmuZj37lR9+VhmUt9MYJG7ry+q4/Vx\naHPfKtp8OvDdonZ2A9uk6ljs7s8V1TGBcE8GXa+eYTKbtCk9TKqeoYhkUz2C4bT4+4Wi60sJz9pK\n5S+Vd4yZTQZ2Apa5e1eSGF8vi+VNAsZWqK9cHcQ8O1Vqs7v/IQZIAMxsF0KP9bdV1jHo0j2/phLD\npB3qGYpIRtUjGI4ButPBK2oHRpXJv6FEXmL+Uunp8sbE96XKSOorV0cuVUbVbY69yVsIwe6r5epw\n904gX6qMwdArGDb1XWeoCTQiklX1mE26Hmgws4aiWZTNQFuZ/M1F15L3bWXS0+WtL/pMqfrK1ZFP\nlVFVm81sd0JvsBk4yN3XlKvDzJoIAbfU9+6lpWV8pSwVtXUWhkEnbzOWlpbxrNzQ2XNt7LhRA1LP\nYNsa2ril6F4U6F4U6F7Urh7BcEn8vQO9hw13pO8wYpJ/h6JrOwJr3X2VmS0BpphZLj4LxMwaCbM2\nn3f3VjNrK1NGUt8SYI8S6cQ8I6tps5m9lRAIlwOHuPvSVN4lwLv6qaNfy5evqZSlomWpMtatbWf5\n8jWsWV3orL7S2jYg9QymlpbxQ76NW4ruRYHuRYHuRUEt/yjYpGFSM5tuZl8ws2vMbBczO8zMplT5\n8ceAtcBBqfJ2A3aj70QZgIXAgUXXZhMmyRB/NwHTU+mzCD2uJM99RfXlYpnJhJiFwH5mlh6unA24\nu6+o0OYF8f2ewG3AP4GZRYEwqWN3M9spdW02sBp4tMT3HnDpCTKFYdLUDjTajk1EMqqmnqGZjQR+\nDLwX6ABGENb4nQO80cxmufs/+yvD3Tea2dXAZWb2CqEX9W3gLnd/MC4zmAy0unsHYX3gOWb2HcLC\n+jmEZQ2HxfKWmtmNhLWHcwkB/lrgOnd/KVb7DeBmM3sUuBM4izCLM1kQfxNwETDfzM4D3kxYjnFK\nFW1+KJZxHWEo9CNAs5lNjdc73f0Vd7/fzB4AfmpmpwPbE54nfj0+Oxx0JSfQ9HpmqNmkIpJNtfYM\nLwLeCfwHYZZm0q34KGHB+sVVlvN5QlC9HrgDeAZ4f0ybQZh4Mh0g7vByOGF5wyPAx4Bj3X1Bqry5\nhN7fLYTAdnvMRyzjd8BJwCeBh4E9gTnu3hrTN8Q6JgAPApcA56Z3uemvzWb2euBthGFPj+1Pfu5P\nlXEU8DKhBzwPuNbdL6zynm22XsGwxHZsmk0qIlmVy+er7w2Y2QvAxe5+dXwu1wHs5+6PmNnRwOXu\nvv0gtTXL8gPxDODRv6/gW794HIDzjtuP1+wwgdXrNvKJby0E4Jg5e/Cvb5vWXxF1p+chBboXBboX\nBboXBS0t43OVcwW19gwnA/8ok7aC0LOSIarXDjSNpbZjU89QRLKp1mD4Z8LzulLeBfxl85ojgyk9\nQaZnBxqddC8iUvPSiouBX8SdX35DWId3gJl9GDiNsIm1DFGdqXWGpZ4ZagKNiGRVTT1Dd7+JEPDe\nCnyPMIHmCuAYwnFGPxvwFsqASfcMkyDY0JCjIZfrky4ikiU1L7p39/mEJQgGbEuYRfrXLXUmn2y6\nXusMU8OjTY05Nnbme/UcRUSyZJN3oIlHIclWpKOzb88QoLGxATq79cxQRDKrYjA0sw7Cs8GquPvI\nyrmkHrq6+550n37d2a2eoYhkUzU9w4spBMNRhIXrfyecTP8iYaj03wmnu180CG2UAdL7PMNCzzDp\nJapnKCJZVTEYuvv5yWszm0fY5eU/k02xo0vM7EeEXVhkiEqCXWNq0kzyHnSEk4hkV63PDP8LeF9R\nIExcB/xy85skg6UnGDb23pRhRFPSM9QwqYhkU62L7tcCryuT9hagdfOaI4MpCXYjGnv/Z29s0DCp\niGRbrT3DGwhDou3AzYTTG6YSNqw+H/jagLZOBlRXT8+wdzDsmUCjnqGIZFStwfAzwC6EI5K+m7qe\ni9e+NEDtkkGQnErRVDRMqgk0IpJ1NQVDd28H3mdmexEO0N2GsEH3nZXOMZT6SxbdN5XpGWoCjYhk\n1SYtunf3PxM27ZatSGdPz7DomWHSM9Q6QxHJqFpPur+tUh53f+emN0cGU/JMsKmhaDaphklFJONq\n7RmOpO9uNOOANxJmmv5iIBolg6OzzASaRk2gEZGMq/WZ4cGlrpvZNsBvgb8NQJtkkCTBcIQm0IiI\n9FLrOsOS3P1V4MvAmQNRngyOpOfXZ2mFdqARkYwbkGCYMnWAy5MBVHECjYZJRSSjap1AM6PE5UZg\nZ+AC4OGBaJQMjp4JNH2GSZNnhuoZikg21TqBZiGlj3PKAUuAT2x2i2TQJEc49V1nqKUVIpJttQbD\nQ0pcywOrgcd12v3QlhzuW3YHmk795xORbKo1GO4K3OLurxQnmNn2ZnaMu399YJomA62ru8IONN15\n8vk8uVyuz2dFRIazWifQ/ADYvUzaPoSDgGWIqjSBBgoBU0QkSyr2DM3sfwmL6iE8G/xVPLWi2FSg\nqv1JzayBEDiPA8YDtwKnuvuyMvn3Ay4H9gWeBy5y9+tT6aOBK4Cj4ne6ETjT3dtSeY4BziNsNP4Y\ncLq7L0qlvxa4CphJOIrqSne/bFPabGYHAHe4+6ii66cA3yYMLSfdr053H9nf/Roo5c4zTA+bdnZ1\n9wmWIiLDXTV/610E3BF/AB5KvU9+bgOuBN5XZb0XAMcCHyZs+D0N+HmpjGa2HSHwLCIEwyuBeWZ2\naCrbtcAM4AjgSOBg4JpUGYcC84BLYxlPALeZ2bYxfUSsYxXwduDTwPlmNrfWNpvZO4BfU/re7h3T\ntk/97FTqew+GznIbdTc09MkjIpIlFXuG7v4A8ACAmTUBF7r705taYQw8ZwCnufud8drRwDNmtn+s\nL+1EYKW7JzNVnzKztwJnA7eb2TTgg8Ah7v5QLO+jwF1m9il3fzHmne/u82L6ycDsWPZXCEF8KnCC\nu68H/mZmewDnEAJvVW02s68QZtQ+Cby5xNd/E6HHuHxT79/mKDdM2tSUGibV8goRyaCaxsPc/YTN\nCYTRPoT9TBekyl0MPEvocRWbCdxTdO1u4ID4egbQBdyXSr83XptpZrmY9+5UfflYZlLfTGBRDITp\nOl5vZi2E3mQ1bT6M0Du9ssT3ANgL+GuZtEHV1d1NPnb6+swmTW3c3aFgKCIZVM0zw43AAe7+kJl1\nUHqdYSLv7s0VipwWf79QdH0pYfF+qfyPlMg7xswmE4YZl7l7V5Lo7l1mtiyWNwkYW6a+/VJ1lEon\nlpEMZfbbZnffF8DMjiv+Ema2I+H8xyPM7ILYpgVA0nsdVOnhz3LrDKFw5qGISJZUs7TiYgpB4GL6\nD4bVGAN0p4NX1A6MKpN/Q4m8xPyl0tPljYnvS5WR1DcGKJ4I006Y5JKUUUubS9mLcO/agQ8A2xH2\nc73DzPaNBycPmvTwZ9/ZpL0n0IiIZE01zwwvSL0+fwDqXA80mFlD0SL9ZqCtTP7i3mbyvq1Merq8\n9UWfKVVfuTryqTJqaXMf7v57M2tx99bkmpm9h/APjSOAm/r7fEvL+GqqKevVNYV/C0yaOLpXeZNf\nWtvzevyE0Ztd12Ab6u3bknQvCnQvCnQvardJJ92b2d6EYb4+zxzd/b6+n+hlSfy9A72HHXek7zBk\nkn+Homs7AmvdfZWZLQGmmFkuPgvEzBqBKcDz7t5qZm1lykjqWwLsUSKdmCdZ+lBtm0tKB8L4/iUz\nW0Hp4eFeli9fU201JbWuLgTDDes39ipvXVshbfkra5nQ3LhZdQ2mlpbxm30vhgvdiwLdiwLdi4Ja\n/lFQ0wQaM3ubmT0NPEqYpPKH1M/C+LuSxwgHAR+UKnc3YDf6TpQhlntg0bXZsX7i7yZgeip9FmGI\nM8lzX1F9uVhmMiFmIbCfmaWHPGcD7u4rNqHNfZjZ6Wb2QgzUybVdgRbC7NNBlR7+bGwovR0b6Jmh\niGRTrT3DK4Bu4HjC4veaHzC5+0Yzuxq4zMxeAZYTFqLf5e4PxmUMk4FWd+8grA88x8y+E+ufAxxN\nmLmJuy81sxsJSyDmEgL8tcB17v5SrPYbwM1m9ihwJ3AWMCGWDWGI8iJgvpmdR1gWcTZwSjVtrvKr\n3xLrmGdmXyY8M7wcuCdZrjGY0hNoRjSVn0Cj2aQikkW1BsO3AUe7+683s97Px7qvB0YAvwVOi2kz\nCAHrEEKgWGZmhwPfIswqXQwc6+4LUuXNJSxnuAXoJO5AkyS6++/M7CTCDjSXxnLmJMOW7r4h1vEd\n4EHCZJpz07vcVGhzRe7+tJnNIUya+SPQAfyKEHQHXe+eYfkJNFpnKCJZlMvnqx8WM7PngI+5+/8O\nXpOkhPzmPgN4eulqLrou7D738fe9mbe8bruetMUvreGCHz4EwGnv3Zu37tGyWXUNJj0PKdC9KNC9\nKNC9KGhpGV/1qQO1bkL5HeBcMxtTMacMKemeYVOfYVItrRCRbNuUI5z2Al40syeAdUXpeXc/bEBa\nJgOqVzDUBBoRkV5qDYZGmEmaGDGAbZFBVO0ONOoZikgW1RQM3b3USfeyFehvBxoNk4pI1unguozo\n6BUMew+TNjbqCCcRybaaeoZm1k35vUnzhIXp/wCuKFqWIHXW1e8waapn2K2eoYhkT609w7OAjcBf\nCIfdngKcD/wppv8IeBr4fjzvT4aIXusM+5x0r56hiGRbrRNo9iecCP/eZB/Q6EIzmw9Mcvf/ijus\nnAX8ZIDaKZupszu1A03xqRUNWnQvItlWa8/w3cA1RYEw8UPgPfH1bcAbNqNdMsA6O9M9w97/2XO5\nXM9QqXqGIpJFtQbD1cCeZdLeQGHdYTOFo5NkCEg/CyyeQAOFAKnZpCKSRbUOk94AXGxm7cAvCRtW\nTyH0CC8EfmBmE4DTgYcGsqGyefpbZwhhIX47CoYikk21BsNzCactfDv+JLqB64BPA0cRni3OGYgG\nysDoNUza0Ldn2NTTM9QwqYhkT62L7juA48zsS4RTJbYjHG57r7s/DWBmtwI7ufuG8iXJlpYMkzY1\n5sjlSgXDcE0TaEQkizbppHt3/yfwz+LrZvY6d//HZrdKBlyyzrDUECmknhl2q2coItlT66L7icDF\nhBPfRxJOk4cwEWcs4flhY+lPSz0lO9CUC4YjNIFGRDKs1tmklwMnERbW54A2woL70UBLTJMhKBn+\nLF5wn0iup58tiohkRa3B8F3AF939PcB3gSXu/gHCaRaPE453kiEomRjT1FD6P3mThklFJMNqDYaT\ngfvi678A+wG4+1rg68CRA9c0GUjJ8Gfxwb6J5IxDTaARkSyqNRiuACbE138HpprZ5Ph+CbDTQDVM\nBlZPz7DsMKmWVohIdtUaDO8APmtmOxNmk7YCx8W0dxOCpQxBPT3DSsOk6hmKSAbVGgy/AEwDfhT3\nJ/0y8HUzexk4G/j+ALdPBkhXzzBp6Z6h9iYVkSyrddH9M2a2B3F/Unf/hpm9CMwEHnT3/zsIbZQB\n0FHtBBr1DEUkgyoGQzP7Qpnr/5Z6+zKwq5md5+4XDlTjZOD09AzLPDMs9AwVDEUke6rpGZ5POMW+\ni8Ii+3LyhA27ZYhJhj+Lj29KJNe7tLRCRDKommB4LzAD+BvhsN6fJPuQytYj6fEVH+yb0DCpiGRZ\nxWDo7rOFtY79AAAd10lEQVTMbCfg/cAHgC+Z2Z8Ixzn91N1fqLVSM2sgbOt2HDAeuBU41d2Xlcm/\nH2H3m32B54GL3P36VPpo4ArCiRlNwI3Ame7elspzDHAesAvwGHC6uy9Kpb8WuIrw/LMVuNLdL9uU\nNpvZAcAd7j6q6HoL4bSPOcBG4AfAZ9190CNQZ4UdaJJ1hppAIyJZVNVsUnd/wd0vd/fpwGuBnwEf\nBBab2UIzO83MptZQ7wXAscCHgVmEGao/L5XRzLYjBJ5FhGB4JTDPzA5NZbuW0Hs9grDw/2DgmlQZ\nhwLzgEtjGU8At5nZtjF9RKxjFfB2wlFU55vZ3FrbbGbvAH5N6Xv7S8L+rbMIQfWEWO6g66ywUXey\nGF+L7kUki2pdWoG7L3b3S9397YRt2P4XOB5YYma3V/p8DDxnAJ9x9zvd/VHgaGCmme1f4iMnAivd\n/RPu/pS7XwX8iLCUAzObRgjMp7j7Q+5+L/BR4ENmtkMs42xgvrvPc3cHTib0/k6M6e8DpgInuPvf\n3P0nwNeAc2pps5l9BVgAPFvie08nBOyPuPuT7n5rLP/0WP6gKhzhVG6YNPQMOxQMRSSDag6GRV4F\nlsWfPOFQ30r2AcYRggYQAiwhgMwqkX8mcE/RtbuBA+LrGYTJPfel0u+N12aaWS7mvTtVXz6WmdQ3\nE1jk7uuL6nh9HNrct8o2H0bonV5Z5nssdvfniuqYQLgngyrZgLvsbNK45CKfh25NohGRjKk5GJrZ\ndmZ2opndBrxIeM62gTDsV81Q6bT4u/hZ41Jg5zL5S+UdE7eC2wlY5u5dSWJ8vSyWN4lwvFR/9ZWr\ng5gn2Wau3za7+77ufmeJ71BNHYMq2YC7/HmGhSCpSTQikjVVLbo3synAewmTaA4EOgnP2I4HfhM3\n6q7WGKA7HbyidmBUmfwbSuQl5i+Vni5vTHxfqoykvjGE4FmcnkuVUUubS+nTTnfvNLN8DWVssq4K\n5xmmr3d25Rk56AO3IiJDRzWL7u8mDDN2AL8DPkLtATBtPdBgZg1FsyibCecjlsrfXHQted9WJj1d\n3vqiz5Sqr1wd+VQZtbS5lD51mFkThXMh+9XSMr7KavrK5/M9E2gmjB9VsqxJE0cXXm8zhonjSt3S\noWFz7sVwo3tRoHtRoHtRu2p6hgcC3YTzCscSZkCeYGal8ubd/bAK5S2Jv3eg97DhjvQdRkzy71B0\nbUdgrbuvMrMlwBQzy8VngZhZI2HW5vPu3mpmbWXKSOpbAuxRIp2YZ2SNbS5lCeE8yHJ19Gv58jVV\nVtNXetizvb2jZFkb1m/sef3ysjVsTL0fSlpaxm/WvRhOdC8KdC8KdC8KavlHQTXPDO8B/kDo2Yyo\n8DOyTBlpjwFrgYOSC2a2G7AbfSfKACwkBOS02YRJMsTfTcD0VPosQo8ryXNfUX25WGYyIWYhsJ+Z\npYcrZwPu7is2oc2lLAR2j2s203WsBh6tsoxNkg6G5bdjK/xR0IxSEcmaahbdHzyQFbr7RjO7GrjM\nzF4BlhMWot/l7g/GZQaTgVZ37yCsDzzHzL5DWFg/h7Cs4bBY3lIzu5Gw9nAuIcBfC1zn7i/Far8B\n3GxmjwJ3AmcRZnHOi+k3ARcB883sPODNhOUYp1TT5iq/9/1m9gDwUzM7Hdge+CrwdXfvrPU+1iK9\nkL6aCTRaaygiWbO5Sys21eeBHwPXE85IfIYwOQfCUomlxJ5e3OHlcMLyhkeAjwHHuvuCVHlzCb2/\nWwiB7faYj1jG74CTgE8CDxNO3Zjj7q0xfUOsYwLwIHAJcG56l5sKba7WUYRNze8hBOJrt8TG5r17\nhuVOuu89gUZEJEty+bz+4tsK5DfnGcCKVev51HfuB+CEd+3JrLfs2CfPo39fwbd+8TgA5x23H6/Z\nYcIm1zeY9DykQPeiQPeiQPeioKVlfKXDJXrUq2coW1CvYdKm/negAehSz1BEMkbBMAOqGibttc5Q\nzwxFJFsUDDMg3dNLTqcopmAoIlmmYJgBvXqGZYZJe2/HpmFSEckWBcMM6BUM1TMUEelDwTAD0j29\nxgpHOEHhuCcRkaxQMMyAaibQpIOkZpOKSNYoGGZANduxjdAwqYhkmIJhBtS6HZsm0IhI1igYZkA1\nPcORqVmmGzuKj20UERneFAwzoJpnhs0jGmnIhUDZtmFQ9w0XERlyFAwzoJph0lwux5hR4RCTde0K\nhiKSLQqGGdBVxTApUAiGGzoGvU0iIkOJgmEGdHZXXmcIMDYGQw2TikjWKBhmQGdnoWc4op9gOGbU\nCADWKRiKSMYoGGZAsqNMLgcNZbZjg0LPUMOkIpI1CoYZkEygKTd5JjGmWcOkIpJNCoYZkAyTVgyG\nqWHSfF4L70UkOxQMMyCZQNPfTFIoDJN25/O0a+G9iGSIgmEGJIvuK/cMm3peaxKNiGSJgmEGFIJh\n/z3DZJgU9NxQRLJFwTADqp5A06tnqBmlIpIdCoYZkOxA09jQ/3/usRomFZGMUjDMgELPUMOkIiKl\nKBhmQM8zw6ZaeoYaJhWR7GiqnGXgmVkDcDFwHDAeuBU41d2Xlcm/H3A5sC/wPHCRu1+fSh8NXAEc\nRfhONwJnuntbKs8xwHnALsBjwOnuviiV/lrgKmAm0Apc6e6X1dJmM3sn8FXAgKeAc9391lT6KcC3\ngTyQdNM63X1kVTduE/UEw352nwEYPbLwx0E9QxHJknr1DC8AjgU+DMwCpgE/L5XRzLYjBJ5FhGB4\nJTDPzA5NZbsWmAEcARwJHAxckyrjUGAecGks4wngNjPbNqaPiHWsAt4OfBo438zmVttmM3sj8Gvg\np8A+wM3Ar8zsDaky9o55tk/97NT/rdp81U6gaWjIMbpZxziJSPZs8Z5hDDxnAKe5+53x2tHAM2a2\nv7s/UPSRE4GV7v6J+P4pM3srcDZwu5lNAz4IHOLuD8XyPgrcZWafcvcXY9757j4vpp8MzI5lfwV4\nHzAVOMHd1wN/M7M9gHMIgbeaNp8B3O/uX4nt/IKZzQQ+DvyfeO1NwB3uvnzz72T1ql1nCGGodH17\np4ZJRSRT6tEz3AcYByxILrj7YuBZQo+r2EzgnqJrdwMHxNczgC7gvlT6vfHaTDPLxbx3p+rLxzKT\n+mYCi2IgTNfxejNrIfQmK7V5VrqOVBnp77QX8NcS33FQVbvOEArLKzRMKiJZUo9gOC3+fqHo+lJg\n5zL5S+UdY2aTCcOMy9y9Z/+w+HpZLG8SMLZCfeXqIOZJhjI3pYydAcxsJ2Ab4Agz+6uZPWdm15vZ\nDiW+84DqqnKYFAqbdWtphYhkST2C4RigOx28onZgVJn8G0rkJeYvlZ4ub0x8X6qMpL5ydeRSZVRq\nc7kykvQ3EibOtAMfAI4H9gDuMLPmEu0fMB01DZPGzbr1zFBEMqQewXA90BBnZ6Y1A21l8hcHi+R9\nW5n0dHnriz5Tqr5ydeRTZVRqc7ky2gDc/fdAi7uf7O6Px2eP7yHMPD2iRPsHTNcmDZPqmaGIZEc9\nllYsib93oPew4o70HWZM8hcPJe4IrHX3VWa2BJhiZrn4LBAzawSmAM+7e6uZtZUpI6lvCaGXVpxO\nzJMsfeivzeXa2ZPf3VvTie7+kpmtoPTwcC8tLeMrZSkrjpIyblxzxXK2mzwWgPUbOjerzsE0VNtV\nD7oXBboXBboXtatHMHwMWAscBMwHMLPdgN3oO1EGYCFhSDFtNmGSDPF3EzCdwiSaWYQhziTPfbG+\nH8f6csCBwHdTdXzIzEa5ezLUORtwd19hZqv7afOCVBkHEdYiJg5J0s3sdOBcYJdkuNXMdgVagCdL\nfO9eli9fUylLWR2dYXS3Y2Nn5XK6Qy9yY2c3S19cyYimxk2udzC0tIzfrHsxnOheFOheFOheFNTy\nj4ItHgzdfaOZXQ1cZmavAMsJC9HvcvcH4zKGyUCru3cQ1geeY2bfISysnwMcDRwWy1tqZjcSlkDM\nJQz9Xgtc5+4vxWq/AdxsZo8CdwJnARNi2QA3ARcB883sPODNhOUYp1TR5odiGVcCi8zsfOAG4Bjg\nXygsq7gl1jHPzL4MbEfYSOCeZLnGYOnsrH4CTfH+pBPHDa1gKCIyGOq16P7zhF7a9cAdwDPA+2Pa\nDMIszOkAcYeXwwnLGx4BPgYc6+4LUuXNJfT+biEEtttjPmIZvwNOAj4JPAzsCcxJhi1jb/BwQoB8\nELiEsHvM9ak6+msz7v4kYQec/wT+RFj8f6S7e0x/mhDIdwb+CPwKeJTw3HDQdHfn6c7XMJt0lHah\nEZHsyeXjX5QypOU3ddijo7OLky8L/274z4N2593Td+s3/xNPv8I3f/YYAJ/98Nt43bSJm1TvYNEQ\nUIHuRYHuRYHuRUFLy/jKswYjbdQ9zHV0Fv6xU8s6Q4B17ZpRKiLZoGA4zHXGCTGgYVIRkXIUDIe5\njo50MKw8YjA2daahdqERkaxQMBzmVqwqbLe6zfhSG/z01rtnqGFSEckGBcNh7uVXC8Fw6uTRFfM3\nNTYwckT4Y6GeoYhkhYLhMLcsBsOGXI5tJ1TuGUJqf1IFQxHJCAXDYW7Zq+sA2G7SqKom0ID2JxWR\n7FEwHOaSYdIp21QeIk2M1TFOIpIxCobDWD6f7xkmnbrNmAq5C8boGCcRyRgFw2FsVdtG2jvCJt21\n9AyTYdJ1GiYVkYxQMBzGlqVnktbUM0yeGapnKCLZoGA4jL3cuq7n9dRanhnGYdING7voSu1gIyIy\nXCkYDmPLVqaWVUysblkF9N6fdH1714C3S0RkqFEwHMaSmaTbTax+WQVoFxoRyR4Fw2FsWRwmrWXy\nDGh/UhHJHgXDYSqfz/PyytqXVUDvnqGCoYhkgYLhMLW6bSPtG2tfVgEaJhWR7FEwHKZq3aA7TcOk\nIpI1CobD1MuvFpZVTNmMYVL1DEUkCxQMh6n0aRXb1bCsAmBkU0PPQcDakk1EskDBcJhKguG2E5tr\nWlYBkMvletYaaphURLJAwXCYSoZJa51Jmkg269aWbCKSBQqGw1D6tIpaZ5ImxmqzbhHJEAXDYWj1\nug429Cyr2LSe4aRxzQD844VVtK7eMGBtExEZihQMh6Flr27aBt1ps982DYCNHd3Mv/3vA9IuEZGh\nqqlyloFnZg3AxcBxwHjgVuBUd19WJv9+wOXAvsDzwEXufn0qfTRwBXAU4TvdCJzp7m2pPMcA5wG7\nAI8Bp7v7olT6a4GrgJlAK3Clu19WS5vN7J3AVwEDngLOdfdbU+ktwLeBOcBG4AfAZ919QI+GePal\nNT2vN3WY9A27bsP+e03lgT+/zCNPLeexf6zgLa/bbqCaKCIypNSrZ3gBcCzwYWAWMA34eamMZrYd\nIfAsIgTDK4F5ZnZoKtu1wAzgCOBI4GDgmlQZhwLzgEtjGU8At5nZtjF9RKxjFfB24NPA+WY2t9o2\nm9kbgV8DPwX2AW4GfmVmb0iV8UtgSvz8ccAJsdwB89SSlfz87n8CMG70CFombVowBPjAIa9jdHMj\nAD/+/VNs7NAJFiIyPG3xYBgDzxnAZ9z9Tnd/FDgamGlm+5f4yInASnf/hLs/5e5XAT8Czo7lTQM+\nCJzi7g+5+73AR4EPmdkOsYyzgfnuPs/dHTiZ0Ps7Maa/D5gKnODuf3P3nwBfA86poc1nAPe7+1di\nO78A3Ad8PJYxnRCwP+LuT8Ye4znA6bH8zbZk2Vqu+PnjdHR205DL8d/vfkPNyyrSJo5r5r0HvhaA\nFas28Jv7niWfzw9EU0VEhpR69Az3AcYBC5IL7r4YeJbQYyo2E7in6NrdwAHx9QygixB4EvfGazPN\nLBfz3p2qLx/LTOqbCSxy9/WpMu4GXh+HNvetos2z0nWkykjXsdjdnytKn0C4J2Xd86fn+ecLq+jo\nLN8ze+7lNXzjZ4+yPi6SP+GIPdlnAIY1D9l3J3adOh6AW+5fzJd+uIiFj7/Yb1sAOru6Wb5yPc+8\nuJpnXlzNsy+tZvFLa1ixaj3t6mGKyBBTj2eG0+LvF4quLwV2LpP/kRJ5x5jZZGAnYJm79/wN6+5d\nZrYsljcJGFumvv1SdZRKJ5axUxVtLldGpfSkjoco49IfPQxAY0OOaVPGseO2Y2mZNIptJ47i1dXt\nPOTLeGF5z+NR/uuQ13HA3juUK64mDQ05jn/XnnzthkdY397F4pfX8P3/91fm3/4UO203lu23HUPL\npNG0re+kdfUGXlm9gdY17axu29hvuc0jGpk4biSTxzczecIothnfzKRxzUwaN5KJ45oZP2YEY0eN\nYMyoJhpyuQH5LiIi5dQjGI4ButPBK2oHSu0bNgYontvfHn+PKpOeLi9ZW1CqjKS+MUDx5J12IJcq\no1Kby7WzbLq7d5pZntLfu4+u7jyLX1rD4tQEmWJHztiVw9+xSzXFVW3X7cdzyYn7c/ejS7n7Ty+w\nqm0jGzZ28c+lq/nn0tWbVGZ7RxfLXl3fsx6ynFwORo9sYlRzI80jGhk7ZiQN+TwjRjQyojFsG9fU\n2EBjQ47GxhwNuRy5hvC7IZcjlwsBPZejJ6jmcjlysez0e3LhP3hPxaTep/KXbmchsd/QXWVcz1WR\ncezYZtra2ivmy4Jx45pZu1b3AnQv0o49cq+q89YjGK4HGsysoWgWZTPQViZ/c9G15H1bmfR0eeuL\nPlOqvnJ15FNlVGpzuTLKpptZE+Gvx1Lfu8f/fG4Oi55cyrMvruGZF1ezbOV6Vq5pJ3l6t/OUcey3\n5xTevucUtp+8aesKK5k4rpn3zHwN756+Kw/7cv66+FVeal3HS6+0sXpdByObGpg8YRTbTmhmmwmj\nenp848eM6PmLvTufp219B6vXbWRV20ZWrt3Iq2s20Lq6nZVr2yn1ODKfD/uj9uyR+sq6vplEREoY\n6sFwSfy9A72HDXek7zBikr94zG9HYK27rzKzJcAUM8vFZ4GYWSNh1ubz7t5qZm1lykjqWwLsUSKd\nmGdkFW0u1850+rv6qaOsqZPH5N594Ov6y7JFHbn9RI6sdyNERAZQPSbQPAasBQ5KLpjZbsBu9J0o\nA7AQOLDo2mzCJBni7yZgeip9FqHHleS5r6i+XCwzmRCzENjPzNLDlbMBd/cVFdqcLuMgejsk9Z0W\nArub2U6p9NnAauDRPt9aRES2mFw9psqb2ZcprLNbTliIvs7d/zUuM5gMtLp7h5lNAf5GWL93BWHB\n+qXAYe6+IJZ3A2FG5lxCgP8+8Ad3nxvTDyOs+zsTuBM4C3gPsGfsOY4C/kwIeucBbwa+R1iucX2l\nNsf0NxHWQn4FuAE4Jtbz1ricAzO7lzD0ejqwPfBD4Cp3v3DAbq6IiNSsXovuPw/8GLgeuAN4Bnh/\nTJtBmGU5HSDu8HI4YXnDI8DHgGOTQBjNJfT+bgFuAm6P+Yhl/A44Cfgk8DCwJzDH3Vtj+oZYxwTg\nQeASwu4x16fq6K/NuPuThB1w/hP4E2Hx/5FJIIyOAl4m9BbnAdcqEIqI1F9deoYiIiJDiTbqFhGR\nzFMwFBGRzKvLqRVSnVpP98gKM7sGaHD3k+rdlnqIk8ouJUwmGw38ETjL3f9c14bVQZydfTlhZnYD\n4f+RT7r7i3VtWB3F/ZL/APyru5eaoT/sxQMS/kyYsJjsYJEHZrn7faU+o57h0Fb16R5ZYWZfIkyG\nyqS4LOhXwOuAfyNMNFsF3GFm29SzbXVyCzCRsKzpQMJa35vr2qI6MrMxhEl+Wf+7fW/CrP/tUz87\nEP7hWJJ6hkNU6qSM09z9znjtaOAZM9vf3R+oawO3MDN7DWEG7l7A4jo3p57eArwDeIO7PwVgZscS\nTmF5N+FEl0wws6nAXwgzv5+L174B3GRmE919VV0bWB/fBJ4Ddq93Q+rsTcBf3H15tR/I+r8ehrJa\nT/cY7mYQ/iffm3APsuo5wpKdp1LXki0CM9UzdPeX3f1DqUA4jXA824NZDIRmdgRhl6szqHoX3GHr\nTcBfa/mAeoZDV62newxr7v5jwjpPzKzOramfuDb2t0WXP07Y7P22Ld+iocHMbiJspNFK2PkpU+Ih\n6P9DmF+wss7NGQreBIwys/sJO4U9CXzW3cueDqSe4dBV6+kekkFm9u+ETSK+XrTBQ9Z8HvgXwraH\nt6cO9s6Ka4Bfufvv692Qeos7iu1OmHR4NuHZ+lJggfXzL2kFw6Gr56SMouvlTveQjDGz4wkTqm5w\n90/XuTl15e5/dvdFwAeBRkIPKRPM7DjCY5Wz46VMD5HGHcUmAbPd/d745+J44GlSO5MVUzAcutKn\ne6SVO91DMsTMPkfYg/dqdz++zs2pCzObYmYfSF9z9/XAPykcyJ0FxxEeq7xsZmsIezkD/NbMrq5f\ns+rH3de6e0fqfZ6w1KLsIyYFw6Gr1tM9JCPM7FPAl4DPu/sn6t2eOtoVuMHM3ppcMLOJgBH+4suK\nY4A3EmYavwU4LF6fC3yhXo2qFzN7q5mtMrN9U9caCL3nJ8t9TnuTDmGVTsrIKjO7C/h7Fhfdm9mb\nCZvN/5DwnCxtjbtn5vTjuObyLsIG+ycDnYRTY3YD9s3SvUiLGxEsAQ7O4qL7eJ7tw8BG4DTCY6VP\nA0cQTipaUepz6hkObf2elJFhWf4X3AcI/9/+N2FSQPonU73EOPT1XsJ5oL8hBMZXCUEgk4EwJbP/\nj8RJh+8CnLABwwOEw95nlQuEoJ6hiIiIeoYiIiIKhiIiknkKhiIiknkKhiIiknkKhiIiknkKhiIi\nknkKhiIiknk6wklkKxNPZLiBcMjvSmCX9D6MIlI7BUORrc8ZwHTgQ8ALCoQim0/BUGTrsy2w1N1/\nUe+GiAwX2o5NZCtiZs8QTmuAsP/kPYSTTU4m7GU7Afg3d19oZgcCFwL7AeuAm4Bz3H1Vqry3AF8n\nDLkuBy4GjgaecfeTzOwgwp6fM939vtTn7gY2uvs74/sc8BnCnqnTCPvoXuru30995i7CfpGLgVOA\nFsKGyh9394dT+fYnnMrxDmADcCtwFmHf0SXAbeljq2LdS4AfZ/1cR9l0mkAjsnX5D8Lmwy8C+wP/\nE6+fB3ycMIT6xxgIfw+sAt5HOPj13cCtyYHR8XSDBcA4wqG4FxCCYc+RSFGpfzEXX7smtuH7wJGE\njbO/Z2anFuX7QGzHxwhBd3vgxiQxHrtzN+GA2mOAU4EDgd/EDZivB46Kp5kn5hDO/fxBiXaKVEXD\npCJbEXd/zMyWA+3u/pCZjYlJV7r7r5J88fivJ93931PX/gT8iRCQbgDOJPyD+HB3XxnzPE0IkFUz\nsz2AjwJnufvl8fLtZtYEXGhm8+Lp48T63pmcKmFmE4Afmtne7v4E8DngZeBd7t4Z87QC3zWz3QnB\n9hzgqPgdAD4CLHL35FBbkZopGIoMD48lL8xsNGGI8ZJ4tlviL4QhyjmEQDITWJgEQgB3/4OZLa6x\n7tnx9/8W1fcbwrFS/0LhQOonio5Xej7+Hht/HwD8KgmEsU13Aq9Pfb/7gQ8TDvYdR+gtn11jm0V6\n0TCpyNYvD6xNvd+G8P/254CO1M9GwvPGHWK+yYTnhMWeL3GtP5MJw5pPFdV3R2zbjqm8xecMdsfP\nJn8XbQssq1DfD4A5ZjaFcL5nA4VeosgmUc9QZPhZTQhClwI/K5G+Jv5eQXhmV2zb1Ovk2WBjUZ5x\nhAktEJ5L5gnP9taXKO+Zyk3usYowsaZHnCDzLuDBeDjrT4BvAv8J/Btwc3pSkMimUM9QZJhx97WE\nZ4N7uPsjyQ/wd+ASwsQbgNuBmWbWExDNzEgNSRICaw7YOZVnG+CNqTzJEOi2RfXtRpjNOpb+pSfj\nLAQOLxpunQn8L2Cp7/cL4FjgYOCHFcoXqUg9Q5GtX67Etc8DvzGzHxB6UqOAc4G9gE/GPJcTJp/c\nbmZfJPT+Liwq73HCsoULzKwtXvsMqWFZd3/CzH4CfN/MLiYE4r2Bi4CH3L3SsGu6vguBewnPH68C\nxhNmuN7p7vem8v2AsORjKfC7CuWLVKSeocjWKV/mNQDufitwOPBaQi9qHtAKHOTuf4l5WoFZhLV/\nPwS+AVxNCH5JOd3Ae4GXCM/lLgfmxzLTPgJ8i7AUIlkX+D3gPf20u8+12KOcTQjePyOsgbyVMCSa\n/n4LCMO917m7FkvLZtOiexHpJS7s/727n1TvtpQT11HeBZi7/6Pe7ZGtn4ZJRWSrYWYHA4cAJwC/\nViCUgaJhUhEplqf0cOZQsB1hs4DnCEOyIgNCw6QiIpJ56hmKiEjmKRiKiEjmKRiKiEjmKRiKiEjm\nKRiKiEjmKRiKiEjm/X+M8/akkbsPqgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d192290>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(0)\n",
    "# plt.subplot(211)\n",
    "# plt.plot(t, np.degrees(resp[:,2]))\n",
    "plt.plot(t,resp[:,0])\n",
    "plt.ylim(20,7)\n",
    "plt.title('Beta Motion')\n",
    "plt.xlabel('time (s)')\n",
    "plt.ylabel('Degrees')\n",
    "\n",
    "plt.figure(1)\n",
    "plt.plot(freq, mag)\n",
    "plt.title('FFT')\n",
    "plt.xlabel('frequency')\n",
    "plt.ylabel('Magnitude')\n",
    "plt.xlim(0, 5)\n",
    "freq_index = np.argmax(mag)\n",
    "print(freq[freq_index])\n",
    "\n",
    "# Should return the zero crossings\n",
    "# zeros = get_zero_crossings(t, np.degrees(resp[:,2]))\n",
    "\n",
    "# Get the peaks of the response\n",
    "# localMaxes, localMax_Times, localMins, localMin_Times = get_local_Extrema(t, np.degrees(resp[:,2]))\n",
    "\n",
    "# Calculate the damping ratio from both the local maxes and the local minimums\n",
    "# zeta_maxes = log_dec(localMaxes[0],localMaxes[-1],len(localMaxes)-1)\n",
    "# zeta_mins = log_dec(localMins[0],localMins[-1],len(localMins)-1)\n",
    "\n",
    "# Average the damping ratio from the positive and negative sides of the plot.\n",
    "# zeta_fromLogDec = np.mean([zeta_maxes, zeta_mins]) \n",
    "\n",
    "# print '\\nThe damping ratio found from the positive peaks is {:.4f}'.format(zeta_maxes)\n",
    "# print 'The damping ratio found from the negative peaks is {:.4f}'.format(zeta_mins)\n",
    "# print 'The damping ratio found from the log dec is {:.4f}\\n'.format(zeta_fromLogDec)\n",
    "resp[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nat_freq_high = np.array([(0)])\n",
    "nat_freq_medium = np.array([(0)])\n",
    "nat_freq_low = np.array([(0)])\n",
    "\n",
    "localMaxes, localMax_Times, localMins, localMin_Times = get_local_Extrema(freq,mag)\n",
    "    \n",
    "first_max = np.argmax(localMaxes)\n",
    "localMaxes[first_max] = 0\n",
    "high_mode = localMax_Times[first_max]\n",
    "nat_freq_high = np.append(nat_freq_high, [(high_mode)], axis=0)\n",
    "\n",
    "second_max = np.argmax(localMaxes)\n",
    "localMaxes[second_max] = 0\n",
    "medium_mode = localMax_Times[second_max]\n",
    "nat_freq_medium = np.append(nat_freq_medium, [(medium_mode)], axis=0)\n",
    "\n",
    "third_max = np.argmax(localMaxes)\n",
    "low_mode = localMax_Times[third_max]\n",
    "nat_freq_low = np.append(nat_freq_low, [(low_mode)], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.          0.48381764] [ 0.          0.91758517] [ 0.         0.4504509]\n"
     ]
    }
   ],
   "source": [
    "print(nat_freq_high,\n",
    "nat_freq_medium,\n",
    "nat_freq_low)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def CRAWLAB_fft(data,time,plotflag):\n",
    "    ''' Function to get the FFT for a response\n",
    "    #\n",
    "    # Inputs:\n",
    "    #   time = time array corresponding to the data\n",
    "    #   data = the response data array (only pass a single dimension/state at at time)\n",
    "    #   plotflag = will plot the FFT if nonzero\n",
    "    #   \n",
    "    # Output:\n",
    "    #   fft_freq = an array of the freqs used in the FFT\n",
    "    #   fft_mag = an array of the amplitude of the FFT at each freq in fft_freq\n",
    "    #\n",
    "    # Created: 03/28/14\n",
    "    #   - Joshua Vaughan\n",
    "    #   - joshua.vaughan@louisiana.edu\n",
    "    #   - http://www.ucs.louisiana.edu/~jev9637\n",
    "    ######################################################################################\n",
    "    '''\n",
    "    \n",
    "    from scipy.fftpack import fft\n",
    "    \n",
    "    # correct for any DC offset\n",
    "    offset = np.mean(data) \n",
    "\n",
    "    # Get the natural frequency\n",
    "    sample_time = time[1] - time[0]\n",
    "    n = len(data)\n",
    "\n",
    "    fft_mag = fft((data - offset)*np.hanning(len(data)))\n",
    "    fft_freq = np.linspace(0.0, 1.0/(2.0*sample_time), n/2)\n",
    "    \n",
    "    # Only return the \"useful\" part of the fft\n",
    "    fft_mag = 2.0/n * np.abs(fft_mag[0:n/2])\n",
    "    \n",
    "    if plotflag:\n",
    "        # Plot the relationshiop\n",
    "        #   Many of these setting could also be made default by the .matplotlibrc file\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",
    "        plt.setp(ax.get_ymajorticklabels(),fontsize=18)\n",
    "        plt.setp(ax.get_xmajorticklabels(),fontsize=18)\n",
    "        ax.spines['right'].set_color('none')\n",
    "        ax.spines['top'].set_color('none')\n",
    "        ax.xaxis.set_ticks_position('bottom')\n",
    "        ax.yaxis.set_ticks_position('left')\n",
    "        ax.grid(True,linestyle=':',color='0.75')\n",
    "        ax.set_axisbelow(True)\n",
    "\n",
    "        plt.xlabel('Frequency (Hz)',fontsize=22,labelpad=8)\n",
    "        plt.ylabel('FFT magnitude',fontsize=22,labelpad=10)\n",
    "    \n",
    "        plt.plot(fft_freq, fft_mag, linewidth=2, linestyle='-')\n",
    "        \n",
    "        # Adjust the page layout filling the page using the new tight_layout command\n",
    "        plt.tight_layout(pad=0.5)\n",
    "        plt.show()\n",
    "    \n",
    "    return fft_freq, fft_mag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# First define all the functions we'll need.\n",
    "def log_dec(peak1,peak2,num_cycles):\n",
    "    '''##########################################################################################\n",
    "    # log_dec.py\n",
    "    #\n",
    "    # Script to compute damping ratio using log dec method\n",
    "    #\n",
    "    # Inputs:\n",
    "    #   peak1 = the amplitude of the first peak\n",
    "    #   peak2 = the amplitude of the Nth peak\n",
    "    #   num_cycles = the number of periods between two peaks\n",
    "    # \n",
    "    # Output:\n",
    "    #   zeta = the damping ratio\n",
    "    #\n",
    "    # NOTE: Plotting is set up for output, not viewing on screen.\n",
    "    #       So, it will likely be ugly on screen. The saved PDFs should look\n",
    "    #       better.\n",
    "    #\n",
    "    # Created: 03/28/14\n",
    "    #   - Joshua Vaughan\n",
    "    #   - joshua.vaughan@louisiana.edu\n",
    "    #   - http://www.ucs.louisiana.edu/~jev9637\n",
    "    #\n",
    "    # Modified:\n",
    "    #   *\n",
    "    #\n",
    "    ######################################################################################\n",
    "    '''\n",
    "    import numpy as np\n",
    "\n",
    "    delta = 1./num_cycles*np.log(peak1 / peak2)\n",
    "\n",
    "    zeta = 1./np.sqrt(1 + (2 * np.pi/delta)**2)\n",
    "    \n",
    "    return zeta\n",
    "    \n",
    "\n",
    "def get_local_Extrema(time,data):\n",
    "    ''' # Function to get the local extrema for a response\n",
    "    #\n",
    "    # Inputs:\n",
    "    #   time = time array corresponding to the data\n",
    "    #   data = the response data array (only pass a single dimension/state at at time)\n",
    "    #\n",
    "    # Output:\n",
    "    #   localMaxes = the amplitude of the local maxes\n",
    "    #   localMax_Times = the times of the local maxes\n",
    "    #\n",
    "    # Created: 03/28/14\n",
    "    #   - Joshua Vaughan\n",
    "    #   - joshua.vaughan@louisiana.edu\n",
    "    #   - http://www.ucs.louisiana.edu/~jev9637\n",
    "    ######################################################################################\n",
    "    '''\n",
    "    from scipy import signal\n",
    "    \n",
    "    # Get local maximums\n",
    "    localMax_indexes = signal.argrelextrema(data, np.greater)\n",
    "    localMaxes = data[localMax_indexes]\n",
    "    localMax_Times = time[localMax_indexes]\n",
    "\n",
    "    # Get local minimums\n",
    "    localMin_indexes = signal.argrelextrema(data, np.less)\n",
    "    localMins = data[localMin_indexes]\n",
    "    localMin_Times = time[localMin_indexes]\n",
    "    \n",
    "    return localMaxes, localMax_Times, localMins, localMin_Times\n",
    "    \n",
    "\n",
    "def get_zero_crossings(time,data):\n",
    "    ''' Function to get the local extrema for a response\n",
    "    #\n",
    "    # Inputs:\n",
    "    #   time = time array corresponding to the data\n",
    "    #   data = the response data array (only pass a single dimension/state at at time)\n",
    "    #   \n",
    "    # Output:\n",
    "    #   zeros = an array of the times of the zero crossings\n",
    "    #\n",
    "    # Created: 03/28/14\n",
    "    #   - Joshua Vaughan\n",
    "    #   - joshua.vaughan@louisiana.edu\n",
    "    #   - http://www.ucs.louisiana.edu/~jev9637\n",
    "    ######################################################################################\n",
    "    '''\n",
    "    \n",
    "    # create an empty zeros array\n",
    "    zeros = []\n",
    "    \n",
    "    for index in range(len(time)-1):\n",
    "        if np.sign(data[index]) != np.sign(data[index + 1]):\n",
    "            zeros.append(time[index])\n",
    "    \n",
    "    return zeros"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}