{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Universidade Federal do Rio Grande do Sul (UFRGS)   \n",
    "Programa de Pós-Graduação em Engenharia Civil (PPGEC)   \n",
    "\n",
    "# PEC00025: Introduction to Vibration Theory\n",
    "\n",
    "\n",
    "### Class 05 - Forced vibration in time domain\n",
    "\n",
    "[1.   Introduction](#section_1)  \n",
    "[2.   Impulse response](#section_2)  \n",
    "[3.   General response to forced vibration](#section_3)  \n",
    "[4.   Short transient loads](#section_4)  \n",
    "[4.1. Impulse response approximation](#section_41)  \n",
    "[4.2. Dynamic factor approach](#section_42)  \n",
    "[5.   Assignment](#section_5)  \n",
    "\n",
    "---\n",
    "_Prof. Marcelo M. Rocha, Dr.techn._ [(ORCID)](https://orcid.org/0000-0001-5640-1020)  \n",
    "_Porto Alegre, RS, Brazil_ \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importing Python modules required for this notebook\n",
    "# (this cell must be executed with \"shift+enter\" before any other Python cell)\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from   MRPy import MRPy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction <a name=\"section_1\"></a> \n",
    "\n",
    "So far we have calculated the system response in damped free vibration, \n",
    "which depends solely on the initial conditions with $F(t) = 0$. \n",
    "\n",
    "<img src=\"images/singleDOF.png\" alt=\"SDOF system\" width=\"240px\"/>\n",
    "\n",
    "We call _forced vibration_ the system response to any force function, $F(t)$, \n",
    "which is the solution of the equation:\n",
    "\n",
    "$$ m\\ddot{u} + c\\dot{u} + ku = F(t)$$\n",
    "\n",
    "In this class we shall present the closed-form solution for any given force function,\n",
    "starting with the very fundamental one: the _unit impulse function_.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Impulse response <a name=\"section_2\"></a> \n",
    "\n",
    "We start by re-writing the equilibrium equation, which after division by the system mass is:\n",
    "\n",
    "$$ \\ddot{u} + 2\\eta \\dot{u} + \\omega_{\\rm n}^2 u = F(t)/m $$\n",
    "\n",
    "Now we take the force function as $F(t) = m\\delta(t)$ (the Dirac's delta function at time origin), \n",
    "then we apply Laplace transform on this equation disregarding the initial conditions \n",
    "($u_0 = 0$ and $v_0 = 0$). It gives:\n",
    "\n",
    "$$ \\mathscr{L}\\left\\{ \\ddot{u} + 2 \\eta \\dot{u} + \\omega_{\\rm n}^2 u \\right\\} = \n",
    "   \\mathscr{L}\\left\\{ \\delta(t) \\right\\} $$\n",
    "   \n",
    "$$  s^2 \\bar{u}(s) + 2 \\eta s \\bar{u}(s) + \\omega_{\\rm n}^2 \\bar{u}(s) = 1 $$\n",
    "\n",
    "what solving for $\\bar{u}(s)$ results in:\n",
    "\n",
    "$$ \\bar{u}(s) =  \n",
    "   \\left[\\frac{\\omega_{\\rm D}}{(s + \\eta)^2 + \\omega_{\\rm D}^2}\\right] \n",
    "   \\left(\\frac{1}{\\omega_{\\rm D}} \\right) $$\n",
    "\n",
    "where we have used the same definition we made for free damped vibration: \n",
    "$\\omega_{\\rm D}^2 = \\omega_{\\rm n}^2 - \\eta^2$.\n",
    "\n",
    "Recognizing the sine function transform and using the translation theorem gives\n",
    "a special response called _unit impulse response_, $g(t)$:\n",
    "\n",
    "$$ g(t) = \\frac{1}{\\omega_{\\rm D}} \\; e^{- \\zeta \\omega_{\\rm n} t} \n",
    "          \\sin \\omega_{\\rm D}t  $$\n",
    "\n",
    "This response corresponds to the free vibration response for $u_0 = 0$ and $v_0 = 1$ (m/s), \n",
    "which means that the unit impulse force $m\\delta(t)$ at time origin causes the same response\n",
    "as a unit initial velocity.\n",
    "\n",
    "From the practical viewpoint, the system response to an impulsive\n",
    "force allow us to know all the system mechanical properties by fitting the theoretical function\n",
    "$g(t)$ to the measured dynamic response.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. General response to forced vibration <a name=\"section_3\"></a> \n",
    "\n",
    "Now let us apply Laplace transform to the general equilibrium equation:\n",
    "\n",
    "$$ \\mathscr{L}\\left\\{ \\ddot{u} + 2 \\eta \\dot{u} + \\omega_{\\rm n}^2 u \\right\\} = \n",
    "   \\mathscr{L}\\left\\{ F(t)/m \\right\\} $$\n",
    "\n",
    "where we define $\\mathscr{L}\\left\\{ F(t)/m \\right\\} = \\bar{f}(s)$.\n",
    "Recalling that the Laplace transform of the unit impulse response, $g(t)$, is:\n",
    "\n",
    "$$ \\bar{g}(s) = \\frac{1}{(s + \\eta)^2 + \\omega_{\\rm D}^2}$$\n",
    "\n",
    "results in:\n",
    "\n",
    "$$ \\bar{u}(s) = \\bar{g}(s)\\bar{f}(s) +  \\bar{u}_0(s) $$\n",
    "\n",
    "where:\n",
    "\n",
    "$$ \\bar{u}_0(s) = \\bar{g}(s) \\left[ u_0 s + \\left( v_0 + 2 \\eta u_0 \\right) \\right] $$\n",
    "\n",
    "is the response to the initial conditions, $u_0$ and $v_0$, while the term\n",
    "$\\bar{g}(s)\\bar{f}(s)$ is the response to the force function, $F(t)/m$.\n",
    "\n",
    "With the _convolution theorem_, the inverse Laplace transform finally leads to:\n",
    "\n",
    "$$ u(t) = u_0(t) + \\frac{1}{m} \\int_0^t {g(t - \\tau) F(\\tau) \\, d\\tau}$$\n",
    "\n",
    "which is the general solution of the dynamic equilibrium equation for a sdof system. \n",
    "It implies that the response, $u(t)$, for any given force function, $F(t)/m$, \n",
    "can be obtained by calculating its convolution with the impulse response, $g(t)$,\n",
    "superposed to the response for initial conditions, $u_0(t)$.\n",
    "\n",
    "The numerical approach for directly solving the convolution integral is not recommended. \n",
    "In the next class we shall be presented to a more efficient numerical technique \n",
    "called _the Duhamel integral_.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Short transient loads <a name=\"section_4\"></a> \n",
    "\n",
    "In this section we present a more practical approach to evaluate the system response to a\n",
    "force function with short duration. To illustrate these approaches, we make use of a numerical\n",
    "example, calculated with ``MRPy``. The numerical method will be explained in details next class. \n",
    "By now we are only interested in the results provided by the Python class.\n",
    "\n",
    "### 4.1. Impulse response approximation\n",
    "\n",
    "Let us start assuming a short transient load consisting in a constant force, \n",
    "$F_0$, applied at $t = 0$ and removed at $t = t_{\\rm d}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "256.0 256\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAv8AAAFECAYAAACwHlEJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0NElEQVR4nO3de3gU9d3+8XuBZJNIEkjIEUIINggKooBGPHAUKlQOQq0CKhRLoQqUpkpBiwQfOVlFK1geBYXSCsRHQPEIsUCAcpAolIOgIJGTCeEUEhIICZnfH/yyZdkAu7DLZjLv13XtJfOd2dnP7mfH697Jd2dthmEYAgAAAFDt1fB3AQAAAACuD8I/AAAAYBGEfwAAAMAiCP8AAACARRD+AQAAAIsg/AMAAAAWQfgHAAAALILwDwAAAFgE4R8AAACwCMI/gCqlQ4cOstls/i4D/9+gQYNks9n0448/+ruUKmHVqlWy2WxKS0tzGvf3+3bu3Lmy2WyaO3eu32oAYA6EfwA+Y7PZPLqZSUUIvPBWu3ZtJSQkqFu3bpoyZYp++uknf5dpKitWrHC8lp988om/y7lmBHIAVVEtfxcAoPoaP368y9iECRMUHh6uUaNGVXqfefPmqbi42MeVeU/r1q314IMPSpKKi4uVm5urdevW6YsvvtCECRP08ssva8SIEX6u8upNnjxZY8aMUf369X3+WO+++66k8x8a33nnHcfragZme98CsC7CPwCfuXhqhHQ+/NepU6fSdZLUsGFD3xblZW3atKn0uXz00Ud68sknNXLkSIWEhOjJJ5+8/sV5QVxcnOLi4nz+OPn5+Vq8eLHuvPNOBQUF6ZNPPtHhw4cVExPj88f2BrO9bwFYF9N+AFQplc2dvnD6xMcff6yUlBSFhISofv36GjdunMrLyyVJ7733nm6//XYFBwerYcOGeuWVVyp9DMMw9O677+qee+5RWFiYQkJC1KZNG8eZZ2/o1auXFi9eLEkaM2aMioqKnNYXFxcrLS1NTZs2VVBQkCIiIvSLX/xC69atc9lXWlqabDabVq1apTlz5qhFixYKDg5WUlKS3njjDcdz+utf/+rYX5MmTfSPf/zDZV/ff/+9Ro8erVatWikyMtKx7ZgxY3Tq1CmX7Sub83/hvPdvvvlGP//5zxUaGqrw8HA99NBDV/X9gPnz5+v06dN6/PHH9cQTT6isrEzz5s2rdNuKmvbu3atXXnlFTZo0UXBwsG6++WYtXLhQklRaWqoXXnhBSUlJCgoK0q233qply5a57Kvi/XbmzBmNHj1aCQkJCgoKUosWLTx6P1z8vh00aJB+/etfS5J+/etfVzq9rVGjRmrUqJFb+6tw/PhxDRs2TDExMQoJCdEdd9yhJUuWXLa2rVu36tFHH1VcXJwCAwOVmJioESNG6NixY24/PwDVB2f+AZjGkiVLtHz5cvXu3Vv33HOPPv30U7300ksyDEN169bViy++qF69eqldu3ZatGiRnn32WcXFxWnAgAGOfRiGoccee0zz589XkyZN1L9/fwUGBiojI0NPPvmkvv3220t+aPBUu3bt1K5dO61evVorVqxQjx49JEklJSXq3LmzNmzYoFatWmnUqFHKy8tTenq6li9frvT0dPXp08dlf6+//rpWrVqlXr16qVOnTlq0aJF+//vfKyQkRP/5z3/0f//3f3rwwQfVqVMnLVy4UE888YSSkpJ07733OvaxePFivfPOO+rYsaM6dOig8vJybdiwQVOnTlVmZqZWr16tgIAAt55fVlaW/vKXv6hDhw4aOnSoNm/erA8//FDbtm3T9u3bFRQU5PZr9c477yggIECPPvqoAgMDNWLECL3zzjt69tlnL3mf1NRUbdy4UT169FDNmjW1cOFC9e/fX3Xr1tWbb76p7du3q3v37jpz5ozmz5+vnj17ateuXUpKSnLZ18MPP6ytW7fq4YcfVmlpqd5//309+eSTOnz4sMaOHev286jQu3dv5efn66OPPlKvXr102223ebyPixUXF6tDhw7atm2b2rZtq/bt2+vAgQN65JFH1LVr10rvs3TpUv3qV79SzZo11bNnTyUkJOjbb7/VjBkztGzZMm3cuFF169a95toAmIgBANeRJCMxMfGS69u3b29c/L+mOXPmGJKMgIAA46uvvnKMFxQUGNHR0UZISIgRGxtr/PDDD451+/fvNwIDA41bb73VaV9vv/22Icl48sknjdLSUsd4SUmJ0aNHD0OSkZWVdcXnsXLlSkOSMXTo0MtuN27cOEOSMW7cOMfYiy++aEgyBgwYYJSXlzvG//Of/xh2u92oW7euUVBQ4BgfP368IcmIiIio9DmGh4cbTZo0MfLy8hzrNm7caEgyevbs6VTPwYMHjZKSEpc6J0yYYEgy/vnPfzqNDxw40JBkZGdnuzx3ScbChQudtn/88ccNScaCBQsu+7pcaPPmzS619uvXz5BkrF271mX7ipqSk5OdnvOGDRsMSUadOnWMe++91zh16pRjXXp6uiHJGDlypNO+Kt5vN998s9NrnpOTY8TFxRm1atVyes0rnvv48eMr3c+FKt63c+bMqfR5JyYmXvJYqGx/Fe+DIUOGOI0vW7bM0Y8LH+vo0aNGWFiY0aBBA2Pfvn1O95k/f74hyRg+fHiljw+g+mLaDwDTGDBggO644w7HcmhoqB588EEVFxfrd7/7nRo3buxYl5CQoHvvvVc7duxQWVmZY3zGjBm64YYbNGPGDNWq9d8/fgYGBmrixImSpAULFnit5vj4eEnS0aNHHWNz585VQECApkyZ4jS149Zbb9WgQYN04sQJffTRRy77GjlyZKXP8eTJk3r++ecVFRXlWHfnnXeqcePG+s9//uO0j/r16yswMNBl38OHD5ckffnll24/t3bt2umRRx5xGhs8eLAkadOmTW7v55133pEkPf74446xJ554wmldZS5+zikpKWrcuLHy8/M1ceJE3XDDDY51ffv2VUBAgMvrceG+QkNDHcuxsbFKTU1VWVmZ5s+f7/Zz8aV58+YpMDBQL774otN4165d1blz50q3Lygo0OTJk12+k9CvXz+1atXKMU0KgHUw7QeAadx+++0uYxVfRq1sWkVcXJzOnTunw4cPq379+iouLta2bdsUHx+vKVOmuGxfWloqSdq1a5fXajYMw2m5oKBAe/fuVbNmzdSgQQOX7Tt06KC33npLW7Zs0WOPPea07mqe/8aNG13qmTNnjubOnavt27fr5MmTju9MSPLo8qStWrVyGat4Tvn5+W7to6SkRO+9957q1KnjmBYlSV26dFFcXJzef/99vfHGG6pdu7bLfS/1euzdu9fl9ahZs6aio6N16NChSuu47777Ljm2ZcsWt56LLxUWFio7O1s333yzYmNjXdbfd999+te//uU0tmHDBsd/9+zZ43KfM2fO6OjRozp69Kjq1avnm8IBVDmEfwCmERYW5jJWcfb+cusqQv2JEydkGIYOHTqkCRMmXPJxLv5y7rXIycmRJMcZ6oKCAkm65FVsKoLdyZMnXdZdzfO/8K8e0vm/HsyYMUMJCQnq2bOn4uLiZLfbJZ2/ElNJSYlbz0uSwsPDL1nPuXPn3NrHkiVLdOLECQ0ZMsRRh3Q+rA8YMECvvPKK0tPTK71a0tW8HhXvhYtFR0e7jFX0qLJeXG8VNVRWp1T5++n48eOSpDfffPOy+y4qKiL8AxZC+AdgGRWBsHXr1srKyrouj7lq1SpJckxXqqjh8OHDlW5fMV5ZeL1WeXl5evPNN3Xrrbdq/fr1CgkJcazLzc297AciX6mY1jNr1izNmjXrktv4+lKpeXl5SkhIcBqr6EVlH3K8oUaNGjp79myl6y7+wFHxfsjLy6t0+8reTxX32bZtm5o3b34tpQKoRgj/ACwjNDRUzZo1086dO5Wfn686der49PEyMzO1Zs0aRUdHq1OnTpLOB7LGjRtrz549OnTokMuPZ2VmZkqqfBrPtdq7d68Mw9D999/vFPwlac2aNV5/vCv58ccf9a9//UsxMTGX/EGvjIwMrV+/Xjt37lSzZs18VsuaNWvUv39/lzHp6ntRs2ZNSZf+K0jdunW1bds2lZWVOX3/pKioSLt373baNiwsTElJSdqzZ49yc3Ndpv5U1r+UlBQtXrxY69evJ/wDcOALvwAsZeTIkSouLtaQIUMqnd6TnZ19Vdepv9jSpUvVt29fSdLUqVOdwvbAgQNVWlqqsWPHOn0nYPv27ZozZ47Cw8PVu3fva67hYomJiZKkdevWOc3zP3jwoMaMGeP1x7uSOXPmyDAMDRs2TLNnz6709sc//lHS5b/46w0TJ05UYWGhY/nw4cOaNm2aatWq5fKhwF0RERGSzr++lWnTpo1KS0v13nvvOcYMw9DYsWMrfW8+/vjjOnv2rF544QWn8eXLl7vM95fO/75AaGionn/+ee3YscNlfXFxseN7AQCsgzP/ACxl6NCh2rBhg/7+97/r3//+t+6//37Fx8fr8OHD2rVrlzZu3Kj58+df8seXLpaVleX4hd8zZ84oJydH//73v/XDDz8oODhYb775pgYNGuR0n9GjR+vTTz/VP/7xD+3cuVOdO3fWkSNHlJ6ertLSUs2bN8/pyjPeEhcXp759+2rRokVq06aNOnfurMOHD+uTTz5Rp06dtHfvXq8/5qWUl5c7frzt4tfnQgMGDNCzzz6rf/zjH5o8ebLbv0HgqcaNG6t58+bq27ev4zr/eXl5mjhxotMVljzRtm1bBQcH6/XXX1dBQYHjex8VH7SGDx+uOXPm6De/+Y0yMjIUFRWlNWvWKD8/Xy1btnS5MtHo0aO1ePFizZo1Szt27FC7du104MABvf/++/rFL36hTz/91Gn7qKgoLViwQA8//LBatmypBx54QE2bNtWZM2e0b98+ZWZm6u6779YXX3xxVc8PgDkR/gFYSsUvBXfv3l2zZs3SJ598olOnTik6OlrJycl65ZVXdP/997u9v6+//lpff/21JCkkJEQRERG65ZZbNGTIED3xxBOOq/FcKCgoSCtWrNDUqVOVnp6u1157TSEhIWrXrp2ee+45px/l8ra5c+eqUaNGWrRokaZPn66GDRsqNTVVf/rTnyq9BKivZGRkaP/+/ercufNlP2hFRkaqR48eWrRokT7++ONKf/zMG95//3298MILWrBggY4cOaLk5GRNmjTpmr5rEBERoQ8++EBpaWmaOXOmTp8+Lem/4b9Fixb64osv9Nxzz+mDDz5Q7dq11b17d/3lL39xuYSqJN1www3KzMzU2LFjtWTJEn3zzTe65ZZblJ6erpMnT7qEf0n6xS9+oc2bN+svf/mLvvzyS2VkZOiGG25QgwYN9Otf/9rlilIAqj+bcfF16AAAsIgOHTooMzPT5ZKsAFBdMecfAAAAsAjCPwAAAGARhH8AAADAIpjzDwAAAFgEZ/4BAAAAiyD8AwAAABZB+AcAAAAsgvAPAAAAWAThHwAAALAIwj8AAABgEVUi/E+ePFl33HGHQkNDFR0drd69e+u7775z2sYwDKWlpSk+Pl7BwcHq0KGDduzY4aeKAQAAAPOpEuE/MzNTTz/9tDZs2KCMjAyVlZWpa9euKioqcmzz8ssva9q0aZoxY4Y2bdqk2NhYdenSRYWFhX6sHAAAADCPKvkjX0eOHFF0dLQyMzPVrl07GYah+Ph4jRo1Sn/6058kSSUlJYqJidHUqVM1dOhQP1cMAAAAVH21/F1AZU6ePClJioiIkCRlZ2crNzdXXbt2dWxjt9vVvn17rVu37pLhv6SkRCUlJY7l8vJyHT9+XJGRkbLZbD58BgAAAMD1YRiGCgsLFR8frxo1Lj+xp8qFf8MwlJqaqnvvvVfNmzeXJOXm5kqSYmJinLaNiYnRvn37LrmvyZMna8KECb4rFgAAAKgiDhw4oAYNGlx2myoX/ocPH66tW7dq7dq1LusuPltvGMZlz+CPHTtWqampjuWTJ0+qYcOG+v777x1/VUDVVlpaqpUrV6pjx44KCAjwdzlwAz0zH3pmPvTMfOiZ+ZipZ4WFhUpKSlJoaOgVt61S4X/EiBFaunSpVq9e7fSpJTY2VtL5vwDExcU5xvPy8lz+GnAhu90uu93uMh4REaHIyEgvVg5fKS0tVUhIiCIjI6v8gYfz6Jn50DPzoWfmQ8/Mx0w9q6jPnWntVeJqP4ZhaPjw4Vq8eLFWrFihpKQkp/VJSUmKjY1VRkaGY+zs2bPKzMzU3Xfffb3LBQAAAEypSpz5f/rppzV//nx99NFHCg0NdczxDw8PV3BwsGw2m0aNGqVJkyYpOTlZycnJmjRpkkJCQtS/f38/Vw8AAACYQ5UI/zNnzpQkdejQwWl8zpw5GjRokCRp9OjROn36tJ566imdOHFCKSkpWr58uVtzmwAAAABUkfDvzk8N2Gw2paWlKS0tzfcFAQAAANVQlZjzDwAAAMD3CP8AAACARRD+AQAAAIsg/AMAAAAWQfgHAAAALILwDwAAAFgE4R8AAACwCMI/AAAAYBGEfwAAAMAiCP8AAACARRD+AQAAAIsg/AMAAAAWQfgHAAAALILwDwAAAFgE4R8AAACwCMI/AAAAYBGEfwAAAMAiCP8AAACARRD+AQAAAIsg/AMAAAAWQfgHAAAALILwDwAAAFgE4R8AAACwCMI/AAAAYBGEfwAAAMAiCP8AAACARRD+AQAAAIsg/AMAAAAWQfgHAAAALILwDwAAAFgE4R8AAACwCMI/AAAAYBGEfwAAAMAiCP8AAACARRD+AQAAAIuoMuF/9erV6tGjh+Lj42Wz2fThhx86rR80aJBsNpvT7a677vJPsQAAAIAJVZnwX1RUpJYtW2rGjBmX3OaBBx5QTk6O4/bZZ59dxwoBAAAAc6vl7wIqdOvWTd26dbvsNna7XbGxsdepIgAAAKB6qTLh3x2rVq1SdHS06tSpo/bt22vixImKjo6+5PYlJSUqKSlxLBcUFEiSSktLVVpa6vN6ce0q+kS/zIOemQ89Mx96Zj70zHzM1DNParQZhmH4sJarYrPZtGTJEvXu3dsxlp6ertq1aysxMVHZ2dkaN26cysrK9PXXX8tut1e6n7S0NE2YMMFlfP78+QoJCfFV+QAAAMB1U1xcrP79++vkyZMKCwu77LamCf8Xy8nJUWJiohYuXKg+ffpUuk1lZ/4TEhKUk5OjyMhIb5cNHygtLVVGRoa6dOmigIAAf5cDN9Az86Fn5kPPzIeemY+ZelZQUKB69eq5Ff5NNe3nQnFxcUpMTNTu3bsvuY3dbq/0rwIBAQFVvolwRs/Mh56ZDz0zH3pmPvTMfMzQM0/qqzJX+/HUsWPHdODAAcXFxfm7FAAAAMAUqsyZ/1OnTmnPnj2O5ezsbG3ZskURERGKiIhQWlqa+vbtq7i4OP3444967rnnVK9ePT300EN+rBoAAAAwjyoT/rOystSxY0fHcmpqqiRp4MCBmjlzprZt26Z58+YpPz9fcXFx6tixo9LT0xUaGuqvkgEAAABTqTLhv0OHDrrcd4+XLVt2HasBAAAAqh/TzvkHAAAA4BnCPwAAAGARhH8AAADAIgj/AAAAgEUQ/gEAAACLIPwDAAAAFkH4BwAAACyC8A8AAABYBOEfAAAAsAjCPwAAAGARhH8AAADAIgj/AAAAgEUQ/gEAAACLIPwDAAAAFkH4BwAAACyC8A8AAABYBOEfAAAAsAjCPwAAAGARhH8AAADAIgj/AAAAgEUQ/gEAAACLIPwDAAAAFlHL0zvs3r1b69atU25urmw2m2JiYnT33XcrOTnZF/UBAAAA8BK3w//Jkyf1xBNP6OOPP1Z4eLiio6NlGIaOHDmigoIC9ejRQ/PmzVNYWJgv6wUAAABwldye9jNixAhlZ2dr/fr1OnHihL777jt9//33OnHihNatW6fs7GyNGDHCl7UCAAAAuAZun/lfunSpli1bppSUFJd1KSkpeuutt/TAAw94tTgAAAAA3uPRF35tNttVrQMAAADgf26H/x49emjIkCHKyspyWZeVlaVhw4apZ8+eXi0OAAAAgPe4Hf6nT5+u+Ph43XnnnYqIiFDTpk3VrFkzRUREKCUlRXFxcXrjjTd8WSsAAACAa+D2nP86dero888/165du7R+/Xrl5uZKkmJjY9W2bVs1bdrUZ0UCAAAAuHYeX+e/adOmBH0AAADAhPiFXwAAAMAiCP8AAACARRD+AQAAAIsg/AMAAAAWcVXhPz8/X7Nnz9bYsWN1/PhxSdI333yjQ4cOXXUhq1evVo8ePRQfHy+bzaYPP/zQab1hGEpLS1N8fLyCg4PVoUMH7dix46ofDwAAALAaj8P/1q1b1aRJE02dOlWvvPKK8vPzJUlLlizR2LFjr7qQoqIitWzZUjNmzKh0/csvv6xp06ZpxowZ2rRpk2JjY9WlSxcVFhZe9WMCAAAAVuJx+E9NTdWgQYO0e/duBQUFOca7deum1atXX3Uh3bp100svvaQ+ffq4rDMMQ6+//rqef/559enTR82bN9ff//53FRcXa/78+Vf9mAAAAICVeHyd/02bNumtt95yGa9fv77jh7+8LTs7W7m5ueratatjzG63q3379lq3bp2GDh1a6f1KSkpUUlLiWC4oKJAklZaWqrS01Gv1HTxxWnmFJVfeEB4rKytTdqH01d6jqlXL47cr/ICemQ89Mx96Zj70zHzM1LMiD2bCePxMgoKCHCH6Qt99952ioqI83Z1bKj5UxMTEOI3HxMRo3759l7zf5MmTNWHCBJfxlStXKiQkxGv1zd5VQ4mhhtf2h4vZtOdfm/xdBDxCz8yHnpkPPTMfemY+5ujZ2dPFbm/rcfjv1auXXnzxRb3//vuSJJvNpv3792vMmDHq27evp7vziM1mc1o2DMNl7EJjx45VamqqY7mgoEAJCQnq2LGjIiMjvVbXxyc2a1r/2y5bC65OaWmpMjIy1KVLFwUEBPi7HLiBnpkPPTMfemY+9Mx8zNSzgoICzXrGvW09Dv+vvPKKunfvrujoaJ0+fVrt27dXbm6u2rZtq4kTJ3q6O7fExsZKOv8XgLi4OMd4Xl6ey18DLmS322W3213GAwICvNtEWw0FBAQQ/n3I6z2Dz9Ez86Fn5kPPzIeemY8ZeuZJfR6H/7CwMK1du1YrVqzQN998o/LycrVq1Ur333+/p7tyW1JSkmJjY5WRkaHbb79dknT27FllZmZq6tSpPntcTxD8AQAAUNVd9bcXOnXqpE6dOnmtkFOnTmnPnj2O5ezsbG3ZskURERFq2LChRo0apUmTJik5OVnJycmaNGmSQkJC1L9/f6/VAAAAAFRnboX/N954w+0djhw58qoKycrKUseOHR3LFXP1Bw4cqLlz52r06NE6ffq0nnrqKZ04cUIpKSlavny5QkNDr+rxAAAAAKtxK/y/9tprTstHjhxRcXGx6tSpI+n8L/6GhIQoOjr6qsN/hw4dZBiXvmKOzWZTWlqa0tLSrmr/AAAAgNW59SNf2dnZjtvEiRN12223aefOnTp+/LiOHz+unTt3qlWrVvqf//kfX9cLAAAA4Cp5/Au/48aN0/Tp03XTTTc5xm666Sa99tpr+vOf/+zV4syDa/wDAACg6vM4/Ofk5FT667jnzp3T4cOHvVIUAAAAAO/zOPx37txZQ4YMUVZWlmOOflZWloYOHerTy30CAAAAuDYeh/93331X9evX15133qmgoCDZ7XalpKQoLi5Os2fP9kWNAAAAALzA4+v8R0VF6bPPPtP333+vXbt2yTAMNWvWTE2aNPFFfQAAAAC85Kp/5KtJkyYEfgAAAMBEPA7/gwcPvuz6d99996qLAQAAAOA7Hof/EydOOC2XlpZq+/btys/PV6dOnbxWGAAAAADv8jj8L1myxGWsvLxcTz31lBo3buyVogAAAAB4n8dX+6l0JzVq6A9/+INee+01b+wOAAAAgA94JfxL0g8//KCysjJv7Q4AAACAl3k87Sc1NdVp2TAM5eTk6NNPP9XAgQO9VhgAAAAA7/I4/G/evNlpuUaNGoqKitKrr756xSsBVV82fxcAAAAAXJHH4X/lypW+qAMAAACAj3k8579Tp07Kz893GS8oKOBSnwAAAEAV5nH4X7Vqlc6ePesyfubMGa1Zs8YrRZmP4e8CAAAAgCtye9rP1q1bHf/+9ttvlZub61g+d+6cvvjiC9WvX9+71QEAAADwGrfD/2233SabzSabzVbp9J7g4GBNnz7dq8UBAAAA8B63w392drYMw1Djxo311VdfKSoqyrEuMDBQ0dHRqlmzpk+KBAAAAHDt3A7/iYmJkqTy8nKfFQMAAADAd9wK/0uXLlW3bt0UEBCgpUuXXnbbnj17eqUwAAAAAN7lVvjv3bu3cnNzFR0drd69e19yO5vNpnPnznmrNgAAAABe5Fb4v3CqD9N+AAAAAHPy+Dr/AAAAAMzJrTP/b7zxhts7HDly5FUXAwAAAMB33Ar/r732mls7s9lshH8AAACginIr/GdnZ/u6DgAAAAA+dk1z/g3DkGEY3qoFAAAAgA9dVfh/55131Lx5cwUFBSkoKEjNmzfX7NmzvV0bAAAAAC9y+xd+K4wbN06vvfaaRowYobZt20qS1q9frz/84Q/68ccf9dJLL3m9SAAAAADXzuPwP3PmTM2aNUv9+vVzjPXs2VO33nqrRowYQfgHAAAAqiiPp/2cO3dObdq0cRlv3bq1ysrKvFIUAAAAAO/zOPw/9thjmjlzpsv422+/rQEDBnilKAAAAADe5/G0H+n8F36XL1+uu+66S5K0YcMGHThwQE888YRSU1Md202bNs07VQIAAAC4Zh6H/+3bt6tVq1aSpB9++EGSFBUVpaioKG3fvt2xnc1m81KJ56WlpWnChAlOYzExMcrNzfXq4wAAAADVlcfhf+XKlb6owy233HKLvvzyS8dyzZo1/VYLAAAAYDZXNe3HX2rVqqXY2Fh/lwEAAACYksfh/8yZM5o+fbpWrlypvLw8lZeXO63/5ptvvFbcxXbv3q34+HjZ7XalpKRo0qRJaty48SW3LykpUUlJiWO5oKBAklRaWqrS0lKv1VVebnh1f/iviteV19c86Jn50DPzoWfmQ8/Mx0w986RGm2EYhic779+/vzIyMvTLX/5SMTExLnP7x48f78nu3Pb555+ruLhYTZo00eHDh/XSSy9p165d2rFjhyIjIyu9T2XfE5Ck+fPnKyQkxGu1zdpVQ0Oall95QwAAAMDLiouL1b9/f508eVJhYWGX3dbj8B8eHq7PPvtM99xzzzUVea2Kiop04403avTo0U5XGLpQZWf+ExISlJOTc8kPDFdj6D83663Hbvfa/vBfpaWlysjIUJcuXRQQEODvcuAGemY+9Mx86Jn50DPzMVPPCgoKVK9ePbfCv8fTfurXr6/Q0NCrLs5bbrjhBrVo0UK7d+++5DZ2u112u91lPCAgwKtNrFHDVuXfFGbn7Z7B9+iZ+dAz86Fn5kPPzMcMPfOkPo9/5OvVV1/Vn/70J+3bt8/Tu3pVSUmJdu7cqbi4OL/WAQAAAJiFx2f+27RpozNnzqhx48YKCQlx+aRx/PhxrxV3oWeeeUY9evRQw4YNlZeXp5deekkFBQUaOHCgTx4PAAAAqG48Dv/9+vXToUOHNGnSpEq/8OsrBw8eVL9+/XT06FFFRUXprrvu0oYNG5SYmHhdHh8AAAAwO4/D/7p167R+/Xq1bNnSF/Vc0sKFC6/r4wEAAADVjcdz/ps2barTp0/7ohYAAAAAPuRx+J8yZYr++Mc/atWqVTp27JgKCgqcbgAAAACqJo+n/TzwwAOSpM6dOzuNG4Yhm82mc+fOeacyAAAAAF7lcfhfuXKlL+oAAAAA4GMeh//27dtfct2WLVuupRYAAAAAPuTxnP+LnTx5Un/729/UqlUrtW7d2hs1AQAAAPCBqw7/K1as0GOPPaa4uDhNnz5d3bt3V1ZWljdrAwAAAOBFHk37OXjwoObOnat3331XRUVF+tWvfqXS0lItWrRIN998s69qBAAAAOAFbp/57969u26++WZ9++23mj59un766SdNnz7dl7UBAAAA8CK3z/wvX75cI0eO1O9+9zslJyf7siYAAAAAPuD2mf81a9aosLBQbdq0UUpKimbMmKEjR474sjYAAAAAXuR2+G/btq1mzZqlnJwcDR06VAsXLlT9+vVVXl6ujIwMFRYW+rJOAAAAANfI46v9hISEaPDgwVq7dq22bdumP/7xj5oyZYqio6PVs2dPX9QIAAAAwAuu6Tr/N910k15++WUdPHhQCxYs8FZNAAAAAHzgmn/kS5Jq1qyp3r17a+nSpd7YHQAAAAAf8Er4BwAAAFD1Ef4BAAAAiyD8AwAAABZB+AcAAAAsgvAPAAAAWAThHwAAALAIwj8AAABgEYR/AAAAwCII/wAAAIBFEP4BAAAAiyD8AwAAABZB+AcAAAAsgvAPAAAAWAThHwAAALAIwj8AAABgEYR/AAAAwCII/wAAAIBFEP69wDD8XQEAAABwZYR/AAAAwCII/wAAAIBFmC78/+1vf1NSUpKCgoLUunVrrVmzxt8lAQAAAKZgqvCfnp6uUaNG6fnnn9fmzZt13333qVu3btq/f79f67LZ/PrwAAAAgFtMFf6nTZumJ598Ur/5zW/UrFkzvf7660pISNDMmTP9XRoAAABQ5Zkm/J89e1Zff/21unbt6jTetWtXrVu3zk9VAQAAAOZRy98FuOvo0aM6d+6cYmJinMZjYmKUm5tb6X1KSkpUUlLiWC4oKJAkjVy4RfaQUK/VFlDTptLSUq/tD/9V8bry+poHPTMfemY+9Mx86Jn5mKlnntRomvBfwXbRBHvDMFzGKkyePFkTJkxwGX8wIk8hIae8Wtdnn/3k1f3BWUZGhr9LgIfomfnQM/OhZ+ZDz8zHDD0rLi52e1vThP969eqpZs2aLmf58/LyXP4aUGHs2LFKTU11LBcUFCghIUEdO3ZUZGSkT+uFd5SWliojI0NdunRRQECAv8uBG+iZ+dAz86Fn5kPPzMdMPauY3eIO04T/wMBAtW7dWhkZGXrooYcc4xkZGerVq1el97Hb7bLb7S7jAQEBVb6JcEbPzIeemQ89Mx96Zj70zHzM0DNP6jNN+Jek1NRUPf7442rTpo3atm2rt99+W/v379ewYcP8XRoAAABQ5Zkq/D/yyCM6duyYXnzxReXk5Kh58+b67LPPlJiY6O/SAAAAgCrPVOFfkp566ik99dRT/i4DAAAAMB3TXOcfAAAAwLUh/AMAAAAWQfgHAAAALILwDwAAAFgE4R8AAACwCMI/AAAAYBGEfwAAAMAiCP8AAACARRD+AQAAAIsg/AMAAAAWQfgHAAAALILwDwAAAFgE4R8AAACwCMI/AAAAYBGEfwAAAMAiCP8AAACARRD+AQAAAIsg/AMAAAAWQfgHAAAALILwDwAAAFgE4R8AAACwCMI/AAAAYBGEfwAAAMAiCP8AAACARRD+AQAAAIsg/AMAAAAWQfgHAAAALILwDwAAAFgE4R8AAACwCMI/AAAAYBGEfwAAAMAiCP8AAACARRD+AQAAAIsg/AMAAAAWQfgHAAAALMI04b9Ro0ay2WxOtzFjxvi7LAAAAMA0avm7AE+8+OKLGjJkiGO5du3afqwGAAAAMBdThf/Q0FDFxsb6uwwAAADAlEwz7UeSpk6dqsjISN12222aOHGizp496++SAAAAANMwzZn/3//+92rVqpXq1q2rr776SmPHjlV2drZmz559yfuUlJSopKTEsXzy5ElJ0vHjx31eL7yjtLRUxcXFOnbsmAICAvxdDtxAz8yHnpkPPTMfemY+ZupZYWGhJMkwjCtvbPjR+PHjDUmXvW3atKnS+37wwQeGJOPo0aPXtH9u3Lhx48aNGzdu3KrD7cCBA1fM3zbDcOcjgm8cPXpUR48evew2jRo1UlBQkMv4oUOH1KBBA23YsEEpKSmV3vfiM//5+flKTEzU/v37FR4efm3F47ooKChQQkKCDhw4oLCwMH+XAzfQM/OhZ+ZDz8yHnpmPmXpmGIYKCwsVHx+vGjUuP6vfr9N+6tWrp3r16l3VfTdv3ixJiouLu+Q2drtddrvdZTw8PLzKNxHOwsLC6JnJ0DPzoWfmQ8/Mh56Zj1l65u6JbVPM+V+/fr02bNigjh07Kjw8XJs2bdIf/vAH9ezZUw0bNvR3eQAAAIApmCL82+12paena8KECSopKVFiYqKGDBmi0aNH+7s0AAAAwDRMEf5btWqlDRs2XPN+7Ha7xo8fX+lUIFRN9Mx86Jn50DPzoWfmQ8/Mp7r2zK9f+AUAAABw/ZjqR74AAAAAXD3CPwAAAGARhH8AAADAIgj/AAAAgEVYJvz/7W9/U1JSkoKCgtS6dWutWbPG3yXhEtLS0mSz2ZxusbGx/i4LF1i9erV69Oih+Ph42Ww2ffjhh07rDcNQWlqa4uPjFRwcrA4dOmjHjh3+KRaSrtyzQYMGuRx3d911l3+KhSRp8uTJuuOOOxQaGqro6Gj17t1b3333ndM2HGtVizs941irWmbOnKlbb73V8UNebdu21eeff+5YXx2PMUuE//T0dI0aNUrPP/+8Nm/erPvuu0/dunXT/v37/V0aLuGWW25RTk6O47Zt2zZ/l4QLFBUVqWXLlpoxY0al619++WVNmzZNM2bM0KZNmxQbG6suXbqosLDwOleKClfqmSQ98MADTsfdZ599dh0rxMUyMzP19NNPa8OGDcrIyFBZWZm6du2qoqIixzYca1WLOz2TONaqkgYNGmjKlCnKyspSVlaWOnXqpF69ejkCfrU8xgwLuPPOO41hw4Y5jTVt2tQYM2aMnyrC5YwfP95o2bKlv8uAmyQZS5YscSyXl5cbsbGxxpQpUxxjZ86cMcLDw43//d//9UOFuNjFPTMMwxg4cKDRq1cvv9QD9+Tl5RmSjMzMTMMwONbM4OKeGQbHmhnUrVvXmD17drU9xqr9mf+zZ8/q66+/VteuXZ3Gu3btqnXr1vmpKlzJ7t27FR8fr6SkJD366KPau3evv0uCm7Kzs5Wbm+t0zNntdrVv355jropbtWqVoqOj1aRJEw0ZMkR5eXn+LgkXOHnypCQpIiJCEseaGVzcswoca1XTuXPntHDhQhUVFalt27bV9hir9uH/6NGjOnfunGJiYpzGY2JilJub66eqcDkpKSmaN2+eli1bplmzZik3N1d33323jh075u/S4IaK44pjzly6deum9957TytWrNCrr76qTZs2qVOnTiopKfF3adD5ecepqam699571bx5c0kca1VdZT2TONaqom3btql27dqy2+0aNmyYlixZoptvvrnaHmO1/F3A9WKz2ZyWDcNwGUPV0K1bN8e/W7RoobZt2+rGG2/U3//+d6WmpvqxMniCY85cHnnkEce/mzdvrjZt2igxMVGffvqp+vTp48fKIEnDhw/X1q1btXbtWpd1HGtV06V6xrFW9dx0003asmWL8vPztWjRIg0cOFCZmZmO9dXtGKv2Z/7r1aunmjVrunxCy8vLc/kkh6rphhtuUIsWLbR7925/lwI3VFyZiWPO3OLi4pSYmMhxVwWMGDFCS5cu1cqVK9WgQQPHOMda1XWpnlWGY83/AgMD9bOf/Uxt2rTR5MmT1bJlS/31r3+ttsdYtQ//gYGBat26tTIyMpzGMzIydPfdd/upKniipKREO3fuVFxcnL9LgRuSkpIUGxvrdMydPXtWmZmZHHMmcuzYMR04cIDjzo8Mw9Dw4cO1ePFirVixQklJSU7rOdaqniv1rDIca1WPYRgqKSmptseYJab9pKam6vHHH1ebNm3Utm1bvf3229q/f7+GDRvm79JQiWeeeUY9evRQw4YNlZeXp5deekkFBQUaOHCgv0vD/3fq1Cnt2bPHsZydna0tW7YoIiJCDRs21KhRozRp0iQlJycrOTlZkyZNUkhIiPr37+/Hqq3tcj2LiIhQWlqa+vbtq7i4OP3444967rnnVK9ePT300EN+rNrann76ac2fP18fffSRQkNDHWcfw8PDFRwcLJvNxrFWxVypZ6dOneJYq2Kee+45devWTQkJCSosLNTChQu1atUqffHFF9X3GPPbdYauszfffNNITEw0AgMDjVatWjlddgtVyyOPPGLExcUZAQEBRnx8vNGnTx9jx44d/i4LF1i5cqUhyeU2cOBAwzDOX4Jw/PjxRmxsrGG324127doZ27Zt82/RFne5nhUXFxtdu3Y1oqKijICAAKNhw4bGwIEDjf379/u7bEurrF+SjDlz5ji24VirWq7UM461qmfw4MGOfBgVFWV07tzZWL58uWN9dTzGbIZhGNfzwwYAAAAA/6j2c/4BAAAAnEf4BwAAACyC8A8AAABYBOEfAAAAsAjCPwAAAGARhH8AAADAIgj/AAAAgEUQ/gEAAACLIPwDABzS0tJ02223+e3xx40bp9/+9rdubfvMM89o5MiRPq4IAKoXfuEXACzCZrNddv3AgQM1Y8YMlZSUKDIy8jpV9V+HDx9WcnKytm7dqkaNGl1x+7y8PN14443aunWrkpKSfF8gAFQDhH8AsIjc3FzHv9PT0/XCCy/ou+++c4wFBwcrPDzcH6VJkiZNmqTMzEwtW7bM7fv07dtXP/vZzzR16lQfVgYA1QfTfgDAImJjYx238PBw2Ww2l7GLp/0MGjRIvXv31qRJkxQTE6M6depowoQJKisr07PPPquIiAg1aNBA7777rtNjHTp0SI888ojq1q2ryMhI9erVSz/++ONl61u4cKF69uzpNPbBBx+oRYsWCg4OVmRkpO6//34VFRU51vfs2VMLFiy45tcGAKyC8A8AuKwVK1bop59+0urVqzVt2jSlpaXpwQcfVN26dbVx40YNGzZMw4YN04EDByRJxcXF6tixo2rXrq3Vq1dr7dq1ql27th544AGdPXu20sc4ceKEtm/frjZt2jjGcnJy1K9fPw0ePFg7d+7UqlWr1KdPH134B+s777xTBw4c0L59+3z7IgBANUH4BwBcVkREhN544w3ddNNNGjx4sG666SYVFxfrueeeU3JyssaOHavAwED9+9//lnT+DH6NGjU0e/ZstWjRQs2aNdOcOXO0f/9+rVq1qtLH2LdvnwzDUHx8vGMsJydHZWVl6tOnjxo1aqQWLVroqaeeUu3atR3b1K9fX5Ku+FcFAMB5tfxdAACgarvllltUo8Z/zxXFxMSoefPmjuWaNWsqMjJSeXl5kqSvv/5ae/bsUWhoqNN+zpw5ox9++KHSxzh9+rQkKSgoyDHWsmVLde7cWS1atNDPf/5zde3aVb/85S9Vt25dxzbBwcGSzv+1AQBwZYR/AMBlBQQEOC3bbLZKx8rLyyVJ5eXlat26td577z2XfUVFRVX6GPXq1ZN0fvpPxTY1a9ZURkaG1q1bp+XLl2v69Ol6/vnntXHjRsfVfY4fP37Z/QIAnDHtBwDgVa1atdLu3bsVHR2tn/3sZ063S11N6MYbb1RYWJi+/fZbp3GbzaZ77rlHEyZM0ObNmxUYGKglS5Y41m/fvl0BAQG65ZZbfPqcAKC6IPwDALxqwIABqlevnnr16qU1a9YoOztbmZmZ+v3vf6+DBw9Wep8aNWro/vvv19q1ax1jGzdu1KRJk5SVlaX9+/dr8eLFOnLkiJo1a+bYZs2aNbrvvvsc038AAJdH+AcAeFVISIhWr16thg0bqk+fPmrWrJkGDx6s06dPKyws7JL3++1vf6uFCxc6pg+FhYVp9erV6t69u5o0aaI///nPevXVV9WtWzfHfRYsWKAhQ4b4/DkBQHXBj3wBAKoEwzB01113adSoUerXr98Vt//000/17LPPauvWrapVi6+wAYA7OPMPAKgSbDab3n77bZWVlbm1fVFRkebMmUPwBwAPcOYfAAAAsAjO/AMAAAAWQfgHAAAALILwDwAAAFgE4R8AAACwCMI/AAAAYBGEfwAAAMAiCP8AAACARRD+AQAAAIsg/AMAAAAW8f8ACPIXsXxRkIYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 900x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "td = 1.00                       # load duration (s)\n",
    "F0 = 10.#/td                       # constant force amplitude (N)\n",
    "Td = 32.                        # total integration time (s)\n",
    "N  = 8192                       # number of time step for Td discretization\n",
    "\n",
    "F  = MRPy.zeros(NX=1, N=N, Td=Td)    # creates a time series filled with zeroes\n",
    "\n",
    "i0 = int(td*F.fs)               # number of time steps to reach td\n",
    "F[0, i0:] = F0                 # fills constant value up to td\n",
    "\n",
    "print(F.fs, i0)\n",
    "\n",
    "F.plot_time(figsize=(9,3), axis_t=[0, Td, -5, 20]);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we use ``MRPy`` to calculate the response of a system with $m = 1$kg, natural vibration frequency\n",
    "$f_{\\rm n} = 1$Hz, and damping ratio of critical $\\zeta = 1$%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "39.47841760435743\n",
      "Maximum response amplitude:    491mm\n",
      "0.15915494309189535\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwwAAAFECAYAAACK67lTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB5FklEQVR4nO3dd3xT5f4H8E/SpmnTvfeEAgXKlqWCTAVluBVlCHLlui5yfyro9Qreq+h1ccVxHQwHywFORKpMZZa9Z+neK23Tpmlzfn/E1Iasc5LUNvB5v168tCcnT5/mSdrzPc/zfb4yQRAEEBERERERWSBv7w4QEREREVHHxYCBiIiIiIisYsBARERERERWMWAgIiIiIiKrGDAQEREREZFVDBiIiIiIiMgqBgxERERERGQVAwYiIiIiIrKKAQMREREREVnFgIGI3N4NN9wAmUzW3t2g382YMQMymQyXLl1q7650CNu2bYNMJsPChQtNjrf3+3blypWQyWRYuXJlu/WBiNwDAwYi6lBkMpmkf+7EeOHY+p+fnx/i4+Mxbtw4vPzyyygoKGjvbrqVLVu2tLyW33//fXt3x2m8iCeijsizvTtARNTa888/b3Zs0aJFCAwMxNy5cy0+55NPPoFGo2njnrlO//79ccsttwAANBoNioqKsGvXLmzatAmLFi3Cf/7zHzz22GPt3EvHLV68GPPnz0dsbGybf6/ly5cDMASay5Yta3ld3YG7vW+J6OrFgIGIOpTLl20AhoAhKCjI4mMAkJCQ0LadcrEBAwZY/Fm++eYbzJo1C48//jhUKhVmzZr153fOBaKjoxEdHd3m36eqqgrr16/HwIED4e3tje+//x7FxcWIjIxs8+/tCu72viWiqxeXJBGR27O0Frz10o7vvvsOgwYNgkqlQmxsLJ577jno9XoAwKpVq9C3b1/4+PggISEBr732msXvIQgCli9fjmuvvRYBAQFQqVQYMGBAyx1uV5g0aRLWr18PAJg/fz7q6upMHtdoNFi4cCG6desGb29vhISE4Oabb8auXbvM2lq4cCFkMhm2bduGFStWID09HT4+PkhOTsZbb73V8jP997//bWmvS5cu+PTTT83aOnv2LJ566in069cPoaGhLefOnz8ftbW1ZudbymFovY7/4MGDuPHGG+Hv74/AwEDceuutDuU7rF69GvX19Zg6dSqmTZuGpqYmfPLJJxbPNfbp4sWLeO2119ClSxf4+Pige/fuWLt2LQBAp9Phn//8J5KTk+Ht7Y1evXrhp59+MmvL+H5raGjAU089hfj4eHh7eyM9PV3S++Hy9+2MGTPwwAMPAAAeeOABi0vvkpKSkJSUJKo9o4qKCsyZMweRkZFQqVS45pprsGHDBpt9O3r0KO655x5ER0fDy8sLiYmJeOyxx1BeXi765yOiKwdnGIjoirZhwwZs3rwZkydPxrXXXosffvgB//73vyEIAoKDg/HCCy9g0qRJGDZsGL766is8+eSTiI6Oxn333dfShiAIuP/++7F69Wp06dIFU6ZMgZeXFzIyMjBr1iycPHnSaqAh1bBhwzBs2DDs2LEDW7ZswYQJEwAAWq0Wo0aNwp49e9CvXz/MnTsXJSUlWLduHTZv3ox169bhtttuM2tvyZIl2LZtGyZNmoSRI0fiq6++wt/+9jeoVCocOXIEX3zxBW655RaMHDkSa9euxbRp05CcnIzrrruupY3169dj2bJlGDFiBG644Qbo9Xrs2bMHr7zyCrZv344dO3ZAoVCI+vkyMzPx6quv4oYbbsBDDz2EQ4cO4euvv8axY8dw/PhxeHt7i36tli1bBoVCgXvuuQdeXl547LHHsGzZMjz55JNWnzNv3jzs3bsXEyZMgIeHB9auXYspU6YgODgY77zzDo4fP47x48ejoaEBq1evxsSJE3H69GkkJyebtXXnnXfi6NGjuPPOO6HT6fD5559j1qxZKC4uxoIFC0T/HEaTJ09GVVUVvvnmG0yaNAl9+vSR3MblNBoNbrjhBhw7dgxDhgzB8OHDkZubi7vvvhtjx461+Jxvv/0Wd911Fzw8PDBx4kTEx8fj5MmTePvtt/HTTz9h7969CA4OdrpvRORGBCKiDg6AkJiYaPXx4cOHC5f/OluxYoUAQFAoFMK+fftajqvVaiEiIkJQqVRCVFSUcOHChZbHcnJyBC8vL6FXr14mbX3wwQcCAGHWrFmCTqdrOa7VaoUJEyYIAITMzEy7P8fWrVsFAMJDDz1k87znnntOACA899xzLcdeeOEFAYBw3333CXq9vuX4kSNHBKVSKQQHBwtqtbrl+PPPPy8AEEJCQiz+jIGBgUKXLl2EkpKSlsf27t0rABAmTpxo0p+8vDxBq9Wa9XPRokUCAOGzzz4zOT59+nQBgJCVlWX2swMQ1q5da3L+1KlTBQDCmjVrbL4urR06dMisr/fee68AQPj111/Nzjf2KTU11eRn3rNnjwBACAoKEq677jqhtra25bF169YJAITHH3/cpC3j+6179+4mr3lhYaEQHR0teHp6mrzmxp/9+eeft9hOa8b37YoVKyz+3ImJiVY/C5baM74PZs+ebXL8p59+ahmP1t+rrKxMCAgIEOLi4oTs7GyT56xevVoAIDz66KMWvz8RXbm4JImIrmj33Xcfrrnmmpav/f39ccstt0Cj0eCvf/0rUlJSWh6Lj4/HddddhxMnTqCpqanl+Ntvvw1fX1+8/fbb8PT8Y2LWy8sLL774IgBgzZo1LutzTEwMAKCsrKzl2MqVK6FQKPDyyy+bLDvp1asXZsyYgcrKSnzzzTdmbT3++OMWf8bq6mo8++yzCA8Pb3ls4MCBSElJwZEjR0zaiI2NhZeXl1nbjz76KADg559/Fv2zDRs2DHfffbfJsZkzZwIA9u/fL7qdZcuWAQCmTp3acmzatGkmj1ly+c88aNAgpKSkoKqqCi+++CJ8fX1bHrv99tuhUCjMXo/Wbfn7+7d8HRUVhXnz5qGpqQmrV68W/bO0pU8++QReXl544YUXTI6PHTsWo0aNsni+Wq3G4sWLzXIs7r33XvTr169lCRcRXT24JImIrmh9+/Y1O2ZMyLW05CM6OhrNzc0oLi5GbGwsNBoNjh07hpiYGLz88stm5+t0OgDA6dOnXdZnQRBMvlar1bh48SLS0tIQFxdndv4NN9yA999/H4cPH8b9999v8pgjP//evXvN+rNixQqsXLkSx48fR3V1dUsOCABJW8H269fP7JjxZ6qqqhLVhlarxapVqxAUFNSyZAsAxowZg+joaHz++ed466234OfnZ/Zca6/HxYsXzV4PDw8PREREID8/32I/rr/+eqvHDh8+LOpnaUs1NTXIyspC9+7dERUVZfb49ddfj19++cXk2J49e1r+e/78ebPnNDQ0oKysDGVlZQgLC2ubjhNRh8OAgYiuaAEBAWbHjLMEth4zBgKVlZUQBAH5+flYtGiR1e9zeYKyMwoLCwGg5U64Wq0GAKu7/xgvBqurq80ec+Tnbz27AhhmKd5++23Ex8dj4sSJiI6OhlKpBGDYwUqr1Yr6uQAgMDDQan+am5tFtbFhwwZUVlZi9uzZLf0ADBf49913H1577TWsW7fO4i5TjrwexvfC5SIiIsyOGcfI0lj82Yx9sNRPwPL7qaKiAgDwzjvv2Gy7rq6OAQPRVYQBAxGRDcaLyP79+yMzM/NP+Z7btm0DgJalVMY+FBcXWzzfeNzSBa+zSkpK8M4776BXr17YvXs3VCpVy2NFRUU2g6i2Ylxy9OGHH+LDDz+0ek5bb0tbUlKC+Ph4k2PGsbAUGLmCXC5HY2OjxccuD1KM74eSkhKL51t6Pxmfc+zYMfTs2dOZrhLRFYQBAxGRDf7+/khLS8OpU6dQVVWFoKCgNv1+27dvx86dOxEREYGRI0cCMFzEpaSk4Pz588jPzzcriLZ9+3YAlpcYOevixYsQBAGjR482CRYAYOfOnS7/fvZcunQJv/zyCyIjI60WacvIyMDu3btx6tQppKWltVlfdu7ciSlTppgdAxwfCw8PDwDWZ1uCg4Nx7NgxNDU1meTT1NXV4dy5cybnBgQEIDk5GefPn0dRUZHZsiRL4zdo0CCsX78eu3fvZsBARC2Y9ExEZMfjjz8OjUaD2bNnW1x6lJWV5VAdgct9++23uP322wEAr7zyiskF+vTp06HT6bBgwQKTHIfjx49jxYoVCAwMxOTJk53uw+USExMBALt27TLJW8jLy8P8+fNd/v3sWbFiBQRBwJw5c/DRRx9Z/Pf3v/8dgO3kZ1d48cUXUVNT0/J1cXEx3njjDXh6epoFEmKFhIQAMLy+lgwYMAA6nQ6rVq1qOSYIAhYsWGDxvTl16lQ0Njbin//8p8nxzZs3m+UvAIb6D/7+/nj22Wdx4sQJs8c1Gk1LngMRXT04w0BEZMdDDz2EPXv24OOPP8Zvv/2G0aNHIyYmBsXFxTh9+jT27t2L1atXWy2odbnMzMyWSs8NDQ0oLCzEb7/9hgsXLsDHxwfvvPMOZsyYYfKcp556Cj/88AM+/fRTnDp1CqNGjUJpaSnWrVsHnU6HTz75xGTHHleJjo7G7bffjq+++goDBgzAqFGjUFxcjO+//x4jR47ExYsXXf49rdHr9S0F+S5/fVq777778OSTT+LTTz/F4sWLRdeIkColJQU9e/bE7bff3lKHoaSkBC+++KLJzlRSDBkyBD4+PliyZAnUanVLHosxOHv00UexYsUKPPjgg8jIyEB4eDh27tyJqqoq9O7d22xHp6eeegrr16/Hhx9+iBMnTmDYsGHIzc3F559/jptvvhk//PCDyfnh4eFYs2YN7rzzTvTu3Rs33XQTunXrhoaGBmRnZ2P79u0YOnQoNm3a5NDPR0TuiQEDEZEdxorR48ePx4cffojvv/8etbW1iIiIQGpqKl577TWMHj1adHsHDhzAgQMHAAAqlQohISHo0aMHZs+ejWnTprXsYtSat7c3tmzZgldeeQXr1q3Dm2++CZVKhWHDhuGZZ54xKbTmaitXrkRSUhK++uorLF26FAkJCZg3bx6efvppi9uttpWMjAzk5ORg1KhRNoOz0NBQTJgwAV999RW+++47iwXtXOHzzz/HP//5T6xZswalpaVITU3FSy+95FTuREhICL788kssXLgQ7733Hurr6wH8ETCkp6dj06ZNeOaZZ/Dll1/Cz88P48ePx6uvvmq2XS0A+Pr6Yvv27ViwYAE2bNiAgwcPokePHli3bh2qq6vNAgYAuPnmm3Ho0CG8+uqr+Pnnn5GRkQFfX1/ExcXhgQceMNuJi4iufDLh8v37iIiIyKobbrgB27dvN9v+lojoSsUcBiIiIiIisooBAxERERERWcWAgYiIiIiIrGIOAxERERERWcUZBiIiIiIisooBAxERERERWcWAgYiIiIiIrGLAQEREREREVjFgICIiIiIiqxgwEBERERGRVQwYiIiIiIjIKgYMRERERERkFQMGIiIiIiKyyrO9O9DR6fV6FBQUwN/fHzKZrL27Q0RERETkEoIgoKamBjExMZDLrc8jMGCwo6CgAPHx8e3dDSIiIiKiNpGbm4u4uDirjzNgsMPf3x8AkJWVhZCQkD/9+wuCgLX7czCyayQiA73/9O/vjnQ6HTZv3oyxY8dCoVC0d3dIBI6Z++GYuR+OmfvhmLkfdxsztVqN+Pj4lutdaxgw2GFchuTv74+AgIA//fsfzKnErpx6ZNcU46Vb0//07++OdDodVCoVAgIC3OLDShwzd8Qxcz8cM/fDMXM/7jpm9pbdM+m5g8s4WYwnb+yK/Mp6CILQ3t0hIiIioqsMA4YO7mxRDbpF+SM+xAd5lfXt3R0iIiIiusowYOjABEGAXhDg6SFH9+hAnCpUt3eXiIiIiOgqw4ChAyup0SIywJDo3DnCD+dLa9u5R0RERER0tWHA0IFdKqtDUpgvgN8DhhLXBgy6Zr1L2yMiIiKiKw8Dhg4sv6oesUE+AIAQXy9Ua3Qua/tAdgWGvrwFR3KrXNYmEREREV15GDB0YHmV9YgL9mmTtlftzcH/7u+Hz/Zkt0n7RERERHRlYMDQgeVX1iO2VcAgl8vQ5IJlRIIgoLy2Ef0TQ1CkbuB2rURERERkFQOGDqy8ToswX2XL15EBSpTUaJ1ut0jdgJggQzJ1TKCPS9okIiIioisTA4YOTBAMswpGMUE+KKhyvhbDkdxq9IoLAgCkRfvjJLdrJSIiIiIrGDB0UJaWCcUG+SDfBQHDxbJadI7wAwCkRQfgdGGN020a1WmbkHmpwmXtEREREVH7YsDQQWkam6FSepocC/dTotQFy4dyyjVIDFEBAJLDfJFdXud0m0YvbjyF//x0Bsfzq13WJhERERG1HwYMHVRZrRZhfl4mx8L8lSitdT5gKKnRItzfkBsR7u+avAgA0DY1o7CqHi9M6oENh/Jd0iYRERERtS8GDB1UWW0jwvyUJsfC/JQoq2l0Sfsymczkv65wLK8afROC0TXS3+VF5oiIiIiofTBg6KAszTAE+ShQXe9cwKBtaoaXh+mwe8plLqn6vP9SJa5JCoFMJoOftydqtU1Ot0lERERE7YsBQwdVXtuIUF/TGQa5XAZnSyYUVDUgJsi0GFxMkA/yK51Ppj5ZqEbP2AAAQM+YQJfmMZTUNOBiKWctiIiIiP5sbhcwvPvuu0hOToa3tzf69++PnTt3inreb7/9Bk9PT/Tp06dtO+giZbVahPkr7Z8oUbG6AVGBpu3GBHmjsLrB6bbrtE3w91YAALpE+uGci5YlCYKAv39+BPO/OoZqjc4lbRIRERGROG4VMKxbtw5z587Fs88+i0OHDuH666/HuHHjkJOTY/N51dXVmDZtGkaNGvUn9dR5lpYkAYBMBuj1jk8zFKsbEOHvbXIsMsAbJTXOBQx6vYDW2RAp4X7IKnXN7kuHcquQHhuI6UOT8O3RApe0SURERETiuFXA8MYbb2DWrFl48MEHkZaWhiVLliA+Ph7vvfeezec99NBDmDJlCoYMGfIn9dR55RaSngEg0McL1fWO32UvrdEiIsC03cgAbxSrnQsYCtWmS53ig32QW6lxqk2jradLcGOPKNzQNRy7zpe5pE0iIiIiEsfT/ikdQ2NjIw4cOID58+ebHB87dix27dpl9XkrVqzAhQsX8Nlnn+Hf//633e+j1Wqh1f6xzahabaiCrNPpoNP9ecthNI06eEAPnc40GTlE5YnCqjr4eTm2u1FhlQbXdwox+VlCVR4orKp36uc7U1iFxBBvkzZ0Tc0uec2O51fhr8OSoJAL0GiboNU2mlTAvpzxe/6Z40XO4Zi5H46Z++GYuR+OmftxtzET20+3CRjKysrQ3NyMyMhIk+ORkZEoKiqy+Jxz585h/vz52LlzJzw9xf2oixcvxqJFi8yOb926FSqVSnrHHVRSIsfGjRvNjxfIsLHiAlIDHVuWdPicHF0aL+B0q5dD2wwcviDHRuGCo93Fb8UyBCqAjZUnWo6Vl8rx3Q8b4eHEzq16ASguliPjp00AAKFGjs++/hFh3naeCCAjI8Nu26svyBHuLeDGOCezyckl7I0ZdTwcM/fDMXM/HDP34y5jptGIWw3iNgGD0eV1AwRBsFhLoLm5GVOmTMGiRYvQpUsX0e0vWLAA8+bNa/larVYjPj4eI0aMQGhoqOMdl+ibikMYP76v2XHNwXz4Kz1xY49IC8+y7+vPDuLWCX3NXrMfPzuI8eP7OdQmAJz++Rxu6hGJ7tEBLcf2NZ9C/2uTzHZlkiKnQoP+QjbGj08DANRk5iFE5YUx3SOsPken0yEjIwNjxoyBQqGwet7mk8UYGViPEwVq9B6ailgn+knOETtm1HFwzNwPx8z9cMzcj7uNmXEljT1uEzCEhYXBw8PDbDahpKTEbNYBAGpqapCZmYlDhw7h0UcfBQDo9XoIggBPT09s3rwZI0eONHueUqmEUmmeO6BQKP60gRcEAXK5zOL3C/P3QVmt1uG+yGVyeHmZJ1PLZXKnfr7imkYkhPqbtBEX4ouSuiYkhjvebk6lFl2iAlra7REbhB1nyzBeRF/tjdmmk6VYNLEHjuVXI+NUGWYPS3G4n+Qaf+bnjFyDY+Z+OGbuh2PmftxlzMT20W2Snr28vNC/f3+zKZ6MjAwMHTrU7PyAgAAcO3YMhw8fbvk3Z84cdO3aFYcPH8agQYP+rK5L1qDTw1vhYfGxYJUClRrXVHu+nOBEkYdqjQ5BKtM3XUyQNwqqnKvvcL6kFp3C/Vq+7hLpj7MlNU61CRh+VnW9DiG+XhiUHIK9WRVOt0lERER0JXKbGQYAmDdvHqZOnYoBAwZgyJAh+OCDD5CTk4M5c+YAMCwnys/PxyeffAK5XI6ePXuaPD8iIgLe3t5mxzua6nodAnwsR3xBKi9UOViLQNPYBB8vy4GIn7cnarRNCPB2PBq+fJlTTJAPMi9VOtweAFworcWkPjEtX/sqPdHQ2OxUm4Z261oCEW+FB/SCAL1esJlMLVZ+VT0eW30Qt/ePw32DEp1uj4iIiKg9uVXAcPfdd6O8vBwvvPACCgsL0bNnT2zcuBGJiYaLssLCQrs1GdyBukGHQKsBgwKVdY7NMBiqR5svRwKAUF8lKmobHQoYrM1MxAT5oLDauboJJTVahF9WwE6A9dwVsQ5mV6J/YnDL18lhvrhYVofOEX42niXOu1vPY/FtvfDSxlOY2DumpZgdERERkTtymyVJRg8//DAuXboErVaLAwcOYNiwYS2PrVy5Etu2bbP63IULF+Lw4cNt30knVdfrrF64B/koHK7DUF2vQ6DKSsDg54VyBwORirpGBFsIRCL9lShysoK0DOYzF0EqhcOzLEZnimvQLdq/5etecYE4nl/tVJsAoG1qRn5VPbpG+eOWXtH45VSJ020SERERtSe3CxiuBtUa6zMMnh5y6B3MNajUNCJYZbndEF8vVDgYMBRWNyAm0HyfU2f6CgC6Zj08LCwRSgzxRXaFc0XhssvrkBjyxza53aMDcKpQ3E4BthzOqcI1SSEAgFFpkdh6xnUBQ622Cav2ZqOmwT32diYiIqIrAwOGDsjWkiRnVGp0CLY2w+DrhfJarcXH7CmpaUBEgIjCCBIVqxsQZSEQSQj1QY6TAUOzXoCnxx9v/4RQFbLLna9MvftiOQanGLbfDfF1rir35RZ9ewIabTMWfXfSZW0SERER2cOAoQMyJD1bTy9x9KZ9labRbCcjI2eWJNnKjZDLZGjWO9bh/Mp6izUcEkJ8kVNe51CbAFDToINKafr6Kj090Nist/IM8Y7nq9ErLrDl66gAbxRWO7dTFAAUVtdD26TH7GEp0DbpUaJ2bqkXERERkVgMGDogdX2T3RkGR7ZArayzPsMQ4qt0eElSeV0jQv3Ma1cAQKif0uGZi4LqeovF1OKCfZBf5fgF88XSOnQON09u9pDLoHMyaGjW66FoNXPRNyEIh3OqnGoTAH4+WYxbekUDACb2jsEPxwqdbtOoWN2AH44WQu9gYEdERERXNgYMHZCtpGfgjy1Qpaqqb7S+JMnP8SVJFXWNCLEywxDur0RJjYMBQ1WDxRmGMD8lyhzsK2DY9jQu2NLMhcqppU6lNVqzwCktOgCnipyvG/Hb+XJc2zkMADC0Uyj2XnRN3YjGJj2eWHcYxwuq8dGvF13SJhEREV1ZGDB0QNX1tnMYglVeqKqTvja+SqNDoJUlSf5KT9Q6EIQAtpckhfsrUergxX1epeUZBg+5zKkic4XVDYgONG83JdwXF0sdX+p0ukiNblH+Jsc6R/jhvJOF5vR6AY3Nevj+vozKV+mJhqZmh5d6tfbj8UKMT4/GUzd2xc5zZWjQOV/jgoiIiK4sDBg6oJoGHfy9recwBKkUqKqXvnyopkGHACvtOlPToLq+0WqAE+6nRKmDMwylNQ1mNRhcoai63nIydYgKuU7MMJwurEFadIDJMZWXJ+qdLDR3sawOKWG+Jse6RQXgdJHzuzptPFaISX1iIJPJMKpbBHacLXW6TaP/bb+A574+Dm0TgxAiIiJ3xoChA9ILpjv4XC5Y5YVKB+sQ2AoMnLhpb7VCcri/4wGDIMDitqqAIUnZ0bvhhdWWd1+KC1Yhv8rxBOWs8jqkhPuaHffylDt10Xw0rwrprRKpASA9NhAnCpwLGBp0zdALaCksN7ZHFDJOFjvVptHPJ4tRpdFhcEoo3tl6wSVtEhERUftgwOCGDIXLHEtQtseRpT62nhLhRMBgS2SAt8NF4TSNzfBTms+0RAd6o8CJgKFE3YAIf/NAJCXcz6mlTkfzqtErLsjkWI+YAJx0MmDYf6kC1yT9Ue06JsjH4XyTy63el4NHRnTC+PQoHMqpdOlSp4M5lbhYWuuy9oiIiMg2BgxuKNDBas/2YgE/b8fzGKxxdIahsUkPTw/rsyHRgd4odDBgsNaqt8IDjU2O75JkbUYkOcwXl8ocDxhyKjRIClWZHHM2QRsAdl8ox9BOYSbHYoK8nZplAYDcCg0i/JXw91ZAJpNhTPdIbD3tmgJ2PxwtxGd7srHwu5M4W+x8MjkRERHZx4ChA7J3YR+k8kKVxCVJzXrB6rIhI0eqPWsam+Dj5WH1cW+Fh0PLccrrtAizslUrAEQFeqNILf3CVq8X4ES6hlW2ZmYSQ1S45GRRuMuXkhnH0pmtUM8W16LrZUna1ySFYH+WczswbTldglFpkS1fj0+PxqYTRU61CRiCyM/2ZOOV23vhP7f3wpsZZ51u00gQBJwoqHZ5wExERHQlYMDQwTQ1662u2zcK8lFIDhjs7bwEGKo9Sw0YbO2QZORIbkRpjdZmwnOUgzMMZXYCEZWDu0WV1TYi1M/y65AY6oucCsdmGDSNTfBRWA7IYoN8UOBgUThBEMxqRgBA34RgHM6tcqhNo9/Ol+HazqEtX4f5KVGp0Tm1sxUAbD9bijHdI6HwkCMq0Bv+3p7IdqKAX2tvZpzF6r05mP1xJqodzA8iIiK6UjFg6GDUDU0IsHNh78guSZWaRgRb2VLVKNiRgMFG0TZn2A0YArxR7EDAUGQl4dkoNsgH+ZXSL8INtR1UFh+L8FeiRO1YbkBOhQYJoZbb7RTuiwsO5kbkVdYj1kItikQnlzo16wXomvVQeZnmiHSN9MMZJ5cQfXekABN6x7R8PbF3LDYec37m4mxxDfIq6/HiremYOzoVb20553SbRhV1jdhwKM/hoohEREQdAQOGDkYtYibA31sBdb20u+BVmkYEWSnaZuTIDENFndZq0TYjby8PyVuLltZoEW4jEIkM8EaRWnrAYKjBYD1gMFSRln7BnFepsVgzAjAsH3L03vqlMvP8BaNOEX64UOJY8u/x/GqkxwaaHZf/XuPC0aVO50pq0CXS3+z4ECeLzema9ajVNpkEkYNSQrA3q9zhNo1W7rqEh0d0/r3NUFwqq4O6wflZhjptEx5ZdRC6JgGPrDqIGhe0SURE1B4YMHQw9qo8A44VLqus01mt8mwUrGqbJUnhDlRmLq3RIiLA+oW9l6ccTc3SL2oNMwyWL+wBIDbYB3mOzDBUWq4ebaTwkDmUy5FdXoeEEPOtWgGgU7gfLji4W9CJAjV6xJgHDIAhodqR1wAADuVUoW9CsNnx3nFBOOLEUqcjuVXoEx9kckzhIUegjwKVTty9r29sRmFVPTpH+LUcm9A7Bj8eK3S4TaO3fjmHR0Z0xl3XxOOxkZ2xdMt5p9sEDHkrS385hwc/3o/MS66p+E1ERGQLA4YOxpBrYL1om6Oq6nV2lySF+HqhQuJ2rWKWJIX5eUmu9lxiZ0mSo+zOMAT5oKBK+syFtSU+RvHBjl2EZ1dokBRmeYYhyomtZS+W1aJTuJ/Fx7rHBOBEQbVD7R7MrkS/hCCz46F+SsnvrdZ+vSwvwmhQcqhTswy/ni/DDV0jTI7d2CMKGSed29WppkGHs8U1uC7VsAvV0M5huFha65LtkD/YeREqpSeW3NMXb205jxwnE+pbO1mgxq4LZU7nmxAR0ZWFAUMHo27Q2c1hcESVphGBIgIGqXdrK+oa7S5JCvNTokzi1qqlNVqEWUkiNpLJZGiWuHTGWpVno+ggH4dqMRSpGxBlY0YkMVTlUIJucXUDIi3UdgCcW+rUoNNb3d2qe3QgThY6VuOhrNb6zJC/t2PbAQOWa1EAhqVOuy84HjBsPVOCEZcFDD5eHvDylDm1hOi7I4W4tV+cybFb+8bhu6POzVyU12qxP6sCM69Ngp/SE4sm9sDrGWecatNo1d5sfPTrRey+UI5nNhxzadDQoGtGU7PjWxYTEVH7YsDQwYjZzchIyh90Q9Kz7QvwQB+F5ArSYpYkhfkpUVYrLRDRNeuh9LS+XauhXS+U10kLRGoamuBvoWibkZ/SE5pG6bsk6fW2q3PHh6gcuhMswHoVbQDw9/aUfBHe1KyH3Mbesp0jHCs0p21qhpen9degd1wgjuVJn7kw5lNcvqMTACSFOr5lrSAIyLWSVH5t5zD8dt65QGRUN9NAZFRaBH52spL2yl2XMOv65JZtdpPDfOHlIXeqzgcAXCitxc6zZXj9zt74+9iuiAn0cTq4Mfp0TzZmf5KJ+5ftxY6zpS5pk4iI/lwMGDoYdb39XZIAwFfpiToJicSVGvs5DJ4ecsnJrtX1jXYDnDB/6TkMYkQEeDu0+9DlNQ0uJ/XGqpjALSFEhVyJS5K0Tc1Q2CheBwApYX7IknixmFdZj/gQ68unfLw8UO9AZeYLJXVWlzkBQM/YQBzNr5Lc7sWyOiSHWc7jkMlk8Pf2dGg24EJprUnuQms3dI3AtjOOLUuqqGuE0lMO38sCU2+FB8L9lQ4vIxMEAYdyqjAkxXRp1n2DE7F6X45DbRot+fkc5o/r1vLZ+MvwFKzZm+NUIUMA2Hq6BKcL1fhk5kB8PHMgVvyW5ZIq3Xq9gDc2n8G05fvw4MeZLttel4iILGPA0MGInWEIUikkrYc27JJkv11HiprZKwgX5uclqdqz2JmTyAAliiXslCS2XblcJmn5RJVGZ3cHKkMytbQ74bkV9UgIsZy/YJQc7ousMmkXYBdKa5Fi48IeMCTW6yQuITlTrDYrBNdaWlQAzhRJ31r1WH6VxR2djHo5OHNxILsS1ySFWHwsNsjHoV24AMNF8qi0CIuPjewWgS0OVr0+V1KL1Eg/s4C3d1wgThWqJS/PMyqpaYAMQFKroEzp6YGbe0Xjh2MFDrUJGALeD3ZcxD9u7g6ZTAalpwdemNQTi3887XCbRm9knEVkoDc+mTkQL0zqgflfHXMq+d2otEaLx9YcwowV+/DIqoMoqXHsPWBJU7PeJTksRETtgQFDByM6YJBYvK2xSQ9vKwXAnCHmGjxM4i5JdY3N8PO2n/gdJXFr1UqNDsF2lk8Z2y2WEODk2dkhCTBcgOkk7uqUXV6HxFDLd9aNUsJ8kSVx+dDF0jp0CrfdbpIDORenC2uQFh1g9fFAlQJqB3IYjuWp0dNGwNA7LgiHHNiB6UB2Jfonmu/oZBTup0SJA0HDrgvluLZTmMXHrk8Nc3hZTsbJYoxpVUHbSCaTYVCy41vM/nisCDf3ijY7fnu/OGw45HjA8NOJYtzcK9okVyY+RIWUcF/svej4cq/j+dUoUjfgvkGJAICYIB88Pa4bXtp4yuE2AcM2uE+sO4xHRnTCygcGYu7oVMxde9glW+z+cLQQ9y/bi+e/PYGZK/c7lCdliV4vYF9WBbaeKUGDA7OCRERiMWDoYNQN9rdVBYBAlZfDCaS2OHJn2R5vhYekpQ0l6gabNRiMIgO8USxhSVJhdb3NHZKMYiQmPudVamzukNSalLyT7HINEq3UYDBKCvPFRYlLki6U1qKznRmGzhF+OC+xxoOtpUNGXp5yyRc2l8rrkGKj3fS4QBzPlz7DUFKjRaSNRPWBySHYmyVt21JBEFBqJ/G7Sa936OJu/6UKXJNseUZkfHo0Njq4Fey2MyUY3iXc7LiPlweSQlU452DBvW8PF2Binxiz47OvT8Hy37IcahMAlv+ahcdHppoc6xMfBKVCjqN5VQ63++pPZ/DXGzqhW5Qh6E2N9Me8MV3w7+9POtwmAGw8VojtZ0vw2axB+O89fbFwQg/8/fMjKHdymWaJugEzP96P386XIau0DjNW7MOhnEqn2jTacbYUT395FP/85rhDny1bNI1NDtd5IaL2w4Chg2lobIa3wv6wSJ1hECtY5YVKkdPmmsYmqzvtOMNelWejiABpd4DtVXk2ignylhQwGKo82w8YQiQWxssur0OSnRkGQ5K2tItPMVvWdo7wx7liaQFDU7PeYmJya6kR/pICEeOFha1lbyovw2sgJRir0jTaDcwHJYdin8SAIadCg0Q7y8gGp4RKDkTKarUI9FFYfX1Twv2QXa6RvCypoKoeIb5Kq7OPt/WLw/pD+ZLaBIDcCg0CfRQWX+MwPyUCfRQOJWoXVtejsVlvMVH9kRGd8d62C5LbBAyzFprGJlzb2XRmaEBSCFReng7vxJVbocG6/bl48db0lk0REkJV+OeE7liw3vGdqGq1TXji88P4x83d8cSYLph5XTI+nDYAS34+5/CWyIDh87bw2xP47XwZnhjTBdOHJmHZr1n4aOdFh9s02nW+DFOX7cWTXxzF9BX78Mqm0y6ZFalvbMYXmbn4z6bT+Hx/Lmq10jetsKZZL0DdoOM2w0RgwNAh2UvKBX7PYah3/XrYYF8vVNaJC0TE7JDkiNJacQFDmK+03Zfs1WAwipFYiyGvsh5xQbYvEgFDLQYpic/5VeL6C0ibuQDsv8c6R/jhvITkVDGVxAEgLTpA0patYmYtgN/HTEIy8aGcKov1IlqLD/FBToW0vJNdF8oxpJN5vYjWhncJx/Yz0pYlbTldgpHdLOdFGA1IDMGBbGl3mH84WohbLCxHMjLubCU1EPkiMxd3DYiz+vi9AxOwZr/0RO2Pd2Vj+tAki49FB/og3F/p0B3x93dcxOOjUi0+Nnd0Kt7bfsGhi8aXN53Gook9zAK9tOgA9EsMxgYHgjEA+Oc3xzF3dBeTpH1/bwWW3N0H//r+pMMXzct/y0JiqAoLxqchKtAbncL98MZdvVFaq8Unuy851CYA/HSiCOsyc/G/+/vjnfv64dNZg9A/IRh/+fSAUzPley+WY+bK/ZDLZJjQOwZennL85ZNMhzcsMCqoqsf/fXEED368H899fRz3L9uLD3dcdHoTgMYmPb45nI+nvzyKJ9YdxvJfs1ySewMYbirsOFuK386XoboNbiYSMWBwU4akZ3G/FBp0zXa3KDUKlXAXXEzRNiOlQvxSFLEzDHKJFa+LqhsQFWB/JkDqkqRCO7UdjBJCpV2A6gXbW7UaSVmaVVnXiCAROTKBPtLyDU4X1dhMeDZKi/bH6ULxS1zsJTwb9Y2XVknakL9geXmPkUwmQ5BKWiXpvRfLMTjFdsDQOcIP50qkLfPZfqYUN3SxHTCM7RGJzSeKJLW783yZ2V311oz5EVJmWpr1Ag7kVGKgleVTgGEJ0ckCtaTq55rGJpwoqMYAG3knD1ybjI93XRLdJgCcK66BSuGBuGDLQX+QygsDk4KRIXFL3F0XyhAT6G2STN7arOuSsf5gvuRE6K1nShDup7SYsB/s64UnRnfBv76TvozqaF4VDuVWYcZlAZlMJsP8m7rhYHalQ7knx/Or8fn+XLx6R2+TncNGd4/Ek2O74rE1hxzayvrnk8X4ePclfDR9AG7vH4e06ABM7huL5TOuwXdHCrHWwZ3Dtp4pwdNfHcWc4SlY8cBA/Peevvhs1iCE+yvx4CeZKKx2LP/kQHYlZqzYh1ptE+aN7YKFE3ogPkSFv607jE/3ZDs8i3GmqAYPrzqAf31/EhdKa3GyQI3564/iyS+OONxXwLDb28s/nsa05fswa+V+PLBiH/762QH8fLLYqRmXWm0TvjqQh6e/PIoHP96P//viCNbuy3E6V0ivF3AopxKr9+bgwx0X8fWhfORKvOFjTZWmEUdyq3AguwL5VfUunXHSNjVL+j3YnlxfUpj+FIE+4nMYKkXukAQYliSJDRgq6rR2i7YZGROfrf1Rbq20Rmv3Lq0jxM4wRPorJSVTNzULNusPGCWEqLDnorgLL12zHh52dp8ySgnzxcWyWlFBy8WyWqTYSXg2kstk0OsFu7tgAYY/WmIChsRQX0nJ1Mfy1JgyKMHueb3jg7D+YB7Gp1u/W97aqUI1/jba8h3l1q5JCsH+SxUY2yPK7rmCIKBCo7P7mZDJZIgL9kFuhQbxdpYvAYaAv0HXbLfwYrcof5wproEgCKJmKbPL6xAb5G33vTuxTwz+t/2C6M/kznOluK5zuM0+yGQy3NgjCptPFGNCb/M8B0u+OpCHO/rH2Ww3OcwXdY1NKKlpQISVgoeX+9/2i3hkRCeb5zxwbTL+8mkmRqdFivo8NOsFvLv1At65r5/VcxQecswdnYo3Ms7ihUk9RfW1vrEZH2y/iGUzBlg9Z1BKKLacLjEkyXc3T5K3pKZBh8UbT+O9+/tZfH1lMhn+fWs6Zn+cibfu9RN1Qwcw5Fm8+MMpvHtfP4vvs/S4QMwZnoL5Xx3Df+/pI+p9CxjyedZl5uKdKebteis88OodvfDs18fh6WF4z4i16XghvjtaiI+mDzC5ySaTyTC5byx6xwdh3roj+Nfknla3ZLbkm8P5+PFYEd6f2h/+rZbpjekeiVHdIrDs1yzM+/wI/j25p9l2zNYIgoBVe3Ow81wpnp/QAzFBpjfDzhbX4Jn1x3BjjyjcfU286NdWEAR8fTgfXx3Ix99Gp+Lpm7q2PLdao8Nne7Px6Z5szB/XzeYmF5fTNeuxak82fjldgkl9YjFvbBdE+BtWCey6UIa5aw+jS6Q/Hh7RSVQep5GmsQmr9uRg29kS9EsIRq+4ICSH+aJIXY+3t5xHSU0DJvWJxfj0aFF/p42qNI1YfzAf286WIshHgU7hfvD0kCG/qh65FRqkRvhjyqAESe8DwPBZ++5IIbaeKUFTs6GAqiAYfs+rvDwxvGs4xvWMMnmfiO3vL6dK8Nv5MlRqGuEhl0EQDO/d7jEBGJYahj7xQaJuQlrDgKEDMfyhF3eulG1Vq0TUYDAK8fVCnsi76+W1jXarMRsZi7eJDRjE/rFXeMihbRI3gyJ2a1lH6lGIER+swhcVeaLOLaiqR2yQuETq5DBfXCytw1ArO/O0dqHUdq2E1uKCfZBXWW9xvfjlThepMbaH/YsTD7kMekEQfVFrL+HZyHDXXtwSqqZmPQRYLgR3uYHJIfgiM1dUwCAmmdxoeJcIbD9bivsHJ9o9d89F+8ucAMMfhm5R/jhVWIPuMfb/kH9/tBC39LJ/sZ4Y6ouCqgbRn7OvDubjuVvS7J43qU8Mnlh3WFTAoNcL+PF4ET6eOdDuufcPSsSqPTl4YkwXu+dml9dBLwh2txn2VXpiRNcI/HCsUFR/vz2Sj5HdIuzueDcgKQRr9+fieH61zZ3AjN7acg4PXp8MlZftP91PjOmCmSv3Y0BisN2d4QRBwPPfnsATY7rYXFbop/TE8xMNuRcfTO1vN3Bq0DXj/748in9N7mmzD0M7heFCSS3e+uW8qCD+dJEab/1yDv+7v7/VC0C5XIYXJ/fE/PVHofCQYVKfWLvtfnM4H1tPl2DJ3X2s/m5IDvPFf+/tgyfWHcaCcWl2x0wQBCz5+RxqGprwzn39LN4EkstlmD0sBbsvlGP2J5lYfFu63d3xahp0eO7r4+gaFYD37rM8Fl0i/bFs+jV4b/sFzPv8CBZN6mH3Qjy3QoMXvj+JfgnBWPHANWavQ6BKgUdGdMad/Rvw0sZTiAzwxuOjUu0GOdvPluKDHRcwuU8sPn5goEl/w/2VmNQnFpP6xOLXc2V4ZNVBDO8SjqlDEm3+vmnQNeOzPdnYeqYE9w1KxKczB5m9Drf2jUOdtgnrD+Vjxop9uC41DFMGJth8nxdU1ePDnReRU67BnQPisWz6AIvvh+P51Vj2axZKaxpwe784jO4eafV9IwgCMrMr8fn+XFRqGjGhdwyW3tvXLHdM3aDDL6eK8cS6I/D39sTEPjG4vnOY1Yt8TWMTfjpRhB+PFUHhIcfIbhF45uY0hLVa+dHUrMfR/GrsOFuKpVvOI9TPC6O6RWJYlzDJQQkDhg6krrHZ7h8Co0AJSc9SZhhCfL1wTOQa4PK6Rru/2IzC/JQoE7lVaaVG3LIZwJD4XFojbuYCEJcfIkV1vU5UoT3A8IuxRORrkF2uQZKIC3UASAn3xW6RSwUulNbi9n7i7rh1jvTHuZIaUQGDYbmXuCAvKtCwHW50oO2ASEzCs5GHXAYZDHd27c3MnC6qQTcRsyEA0DlcfCCyW0T+gtHQzqGY/9VRUQHDz6eKMeu6FFHtju0Rhc0ni0QFDHsulmPOcNt31o1GdA3HtjOluNFO4FReq4VeL4gK+P29FQjzUyJLRJ7KltMluD41XFSQN6RTKN7ddgENuma7W0n/b/tFPDRc3Gt7/+BEzFy5H+N6Rtm8S9fYpMfn+/OwcuY1otp96qaueOrLo1g+/Rqb7/UzRTXIrdDg6Zu62W3TW+GB+eO6YdF3J7Dknr42z/3qYD4SQ3xtLiEz6hYVgDHdI/DutvN4dKT1i3tBELBg/THMvj5Z1B3YqUOS8M9vjuPbIwWYaCMgy63QYNG3J/HOff3sXqTK5TIsvq0XnvziCLw85BhnY/bxqwN52H2xHK/f1cfu748If2+8M6UfHl97GHOGp1i9UdOga8aC9cfQLyFIVPA6pFMoEkJVmP/VUcy6Lhk3dLW8BPFoXhUWbzyNv4/tggFW6sgYyeUyPDKiMw5kV+ChTw7goeEpFtttbNJj2a9ZOJBdgX/c3N3qMjqjiABvLLmnL7adKcGDH2di2pBE3Ngjyuz9e6pQjSU/n0VKuB/enzoAfnbG7LrUMAztFIrvjhZgxvL9uLVvLG7uFW0y1sY7/z+fKsY9AxMsBgqt+So9MXVwIu4bmIBtZ0vwf18cRXSgN6YNSUTnCENdm5oGHfZfqsC3hwvQpBfw4PUp6BMfZLOvPWMDsfi2dKgbdFh/IA8zVuxDemwQBqeEICbAC6X1wG8XynEwpxoHcirRNz4Yj41Mtfn3NMBbgVv7xuHWvnEorK7HN4cLsPzXLMQE+iA9LhBhfl5o0gvIrajH4dxKNOsFjO0ehdfu6m01GPT0kKNfQjD6JRiWcharG/DLqRI8+cVRNAsCOoX7IVwpbkkUA4YORF2vQ4CPuCFReMhFJyKKKSxmJGUnn4q6RglLkrxE12IQBHEXicAf6/ftBQxS1xz6Kj1Rq22y+wsur1IjaockQFqwIqYGg1F8iEp0bsSlsjq7W7UapUb44VBOFUZZ2Pu/NeNsgdifLy06ACcL1HYDBrEJz0adIw25AcZtMa05mGO7/kJrcrkMKi8P1Gmb7F6g7M2qwIu3potqN8BbgQadHo1NepvT5IIgILtcI/p16JcQjKVbzts971xxDZLDfEUve7u5Vwz+9f1JuwHDhkP5uLWv/bu5RvcMTMDa/TlYMM72jMTa/Tl4/c4+otqUyWSY1CcG6w/m21zOVlBVD3WDzu77xchb4YFxPaPwzeEC3G5jmcva/Tm4o3+c6LyxCH9vXJ8ajvWH8q0un9HrBSz+8RQW3ybu/QUAveKCkBCiwg9HCy3W2QAMNxA2HivEh9OsL3G63F0D4vF/XxzFHhv5Om/+fA59E4Jwfar5dr3W/POW7vjb2sMI9fWymFdTWqPFU18exWt39Rb9d8dDLsMrd/TCvM+PAIBZ0CAIAlbuuoQzRTV45fZeoj8PQSov/O/+fnhm/TGcLFDjgWuTTZ57rrgGi747iTnDO+G6VPszv0axQT74cNoALPruBLaeLsEjIzoj2MfwPiqv1WL57vMoVjfgvfv7if6bDgD9E0OwbMYA/Pfnc1izLweT+sQiLToA9Y3N2HWhDL+cKsE9A+MxZ/gASX+nbugagcEpoVjx2yXcv2wvukUFIC7YB1WaRhzNr0ZskA/+OaGH6NlywPA7d1KfWIzrGY2vD+dj7rrDaPp9iW6TXoCvlydu6RWNT2YOlLS8Ri6XYWS3SIzsFomTBWqs2puD3AoNZDJAqfBAv4RgLBifZnOrbUsCvBWYcW0ypg9NwvF8NfZdqsCOMyW4UCxHc1ENRqZF4m+ju4h+bxlFB/pgzvBOmDO8E4qqG3A8vxoVmkbIZTL0iQ/CtCGJopevtRYZ4I0pgxIwZVAC6hubkVepQVGZuKXSDBg6ELFF26Sq1DSKvvsbLCXpWcqSJH+l6Cq/UiYBIvzFba2qbmgSPRMAANFB3iisqkdqpO070fmV4pcOAYDn73Uu7N0pvVSusZmM2pqU4FHXLIi+kOkS6Y/PM3PtnldQ3YCYIPG/ZHvEBGDX+XK7gYjYhGejPnGGxGe7AUN2pailOEb9EoJxMKfS5sWPXi+gVtsk6fPbPzEYB7Irbc5KnCxUo7uEdcIechlig7zt5kd8d1Tc0hqjcH8lGnTNqGnQWZ3GFgQBW8+U4OMH7C8bMuodF4jXN5+xGTidKKhGTJCP3RyO1ib3jcX05ftw14A4qxcV/9t+AX8VOcNidNc18ZixfD8m9omx+BmuadDhpxNF+GTmIEntThuSiBkr9mFM90iL76F1mbkYlhpuN8i+3KMjUzFz5X4MTA4xyzuo0zbhua+P48277d9Vb00mk+GFST3wl08zEebnhc4Rpr8jP92TDa2uGdOGJEnqq6eHHK/d2RuPrTmEKo3OJMg5XaTGom9P4l+Te0r6fQsYfj++fmdvLP7xFH67UIaHhnVCXLAPzhTXYOmW8+gZY7hTLHX2WeXliTfv7oPV+3Jw/0d7MTglFAE+ni2bL7x6Zy/J4wUYAtPFt/XCgewKvLTxFKrrG1FcLMfm2jO4b3CiqKWn1vq7YHwaymu12HyyGJ/uzoa3Qo7+icH4dJa0i+/L+/vXGzrhoWEpuFhWh6LqBgT4eOKRkZ1F/62xxMtTjrsGxOOuAfEQBEHSjUR7uscEYOHEHi5py0gmkyE9LhDpcYHQ6XTYuPECxl+bBIXC+Wu6qEBvUTmKUvl4eSA10h+RPuKuHxgwdCDqenFF26Sq0oi/i+br5SF6x4pqCf0NF1nt2bAMRfwvhahAb1H7+outwWAUG+SDfBEBQ15lPbrYOae1mCAfFFY12F3mk1uhEb3MCvij4J6tQERKIjVgmG0Ss0PQ2aIadJXwGnSLCsDyXy/ZPe9Ynhr3Dbaf8GzUOz4IS7ecx912VoJUikhMbu2apBD8cqrYZsBwqkgt+jNmdEPXcHx7pMBmwPDzyRK7gdXlxnaPwk8nivDg9ZaX2giCgIPZlZhrZRtRa27sEYWfThRbvQt+MKcSveKkJdXJZDKM6R6JjJPFVu+CL7NQqM0ehYccN/aIwg/HCi2uXy+qbkBpjVZU3kBrSk8PTO4bgzX7cixeEL+++SweGdFZ8t1EhYccj41Mxeubz5glQJ8vqUHGyWJJswBGXp5yPHtzGp768gjentIPXr8PTYOuGfO+NGzNKvWOKmCYgV1yd188se4w7h+cgLHdo1Db2IS3fj4HAHhmvP0cFkt8vDzw3v398EbGWaw/mIf0uEDklGvQpBfw33v7iM5tu5yXpxzPT+iBgzmVeHvLeZTVapEQqsK8MV1E53RZIpPJcN+gRNw1IB5H86pRq23CzenRVgs3StE/MQT9E0N+v/jciPHje7nk4jPUT4l7B4r/vSqWXC5D5wg/yUnAYhhmsF3eLEnEbVU7kLaaYajSNCJY5N05aXdZxO2gA/yR9GxPpYS+AuK3FBVb5dkoJlBcLYa8ynrRVZ4B8cuHmvTidl4ySghR2d1CTkpehJHs952SbDlTXCMpaPJVeooKSrPKapEsclkWAEQHetvdRrBY3YAIkTu8GPWICcCJAtu1I3ZfKMdQiTt7dY8OwCk7W8wezKm0Wy/ickM7h2KXjUJjx/PV6B4TIPlu3dgekfjJxratq/fmYooDFyKT+sTiq4OWNwPIrTAUo7O3ptqSu6+Jx7r9uRa3LHx98xk8MqKz5DYB4I7+8dh6usRst6/dF8pR39js8B3gwSmh8FV6YvmvWS1LKHMrNHju6xN4+bZ0yUGIUVp0AB4e0RmzP8nEhkMF2Fcqw0OfHcIDQ5NE5S1YE+6vxEfTB+BCaR0e/CQTT395FNemhuEft3R36k6wwkOOp2/qhiX39MENXSMwf1w3vHVvX4eDhdb6JQTjlTt6YdmMa/D8hB5OBQutKTwMd+qHdwl3SbBA1BExYOhApCTQAoapQDG1DSol5DC0FR8vcX0tq21EmIQLukh/b1FLkqQk5QLiazEUVNVLWo6TEKJCbqXtC3u9XoDUv7cp4Yadkmy5UFor+Q9kXLBhpsWWs0XSAgbAkPBqa1vgZr0AuUwm6cJDJpPBR+GBehuVrw9kV2JAkrj8BSNPDzlkMpnNvbL3X6qQ3K5MJkOEvxLFVt6/+VX1iPBXSl4qoPT0gJ/S0+qM3vpDeZLyDIz8vRVQeXmgyEKBvIq6RtRpm0RtE3u5QB8Fukb5Y9f5MrPH3t12AbOtzJTY463wwP2DE/HuVtPqz7+eK4PKy0Py7IKRh1yGFyb1xPyvjrUEqEfzqvDO1vN4fmJ3h9o0eurGrqjVNuHBjzPx1JdH8Py3J/Dqnb2cvgi9JikE/5vaHx5yGeQAlt7TG0NFLnm0xVvhgUdGdMbyGdfgvfv7Y4SVRF1H+Hsr0Cc+iBfgRB0EA4YORN0gbQ10oMjibVJnLgRBepKwq5TWaBEushgcAAT4eIqqR2GowSB+JiBWZMCga9ZLWqcppnpwkVra8inAsN1fVpmIgEHidHGqiCJjlZpGu1s3Xq57jCHx2ZrzJbUOTW33jA3E8QLru3wZCrZJu7AH/qh4bImuWQ9dsyB6h7PWhncJx/azlqs+bz5RZDfJ2JrR3SPxs4VCY7pmPc6X1EraP721aUOS8NHOi2bHP951CfeKqJdhzUPDUvDutgvQNf9RSfdYXjWamvUOX9gDwLieUSioqse3RwoAGLZC/N/2C3h6nP2dhmyJD1HhxVt74vlvTmDGin1Y+dslvDOln0PvgdZkMhkeH5WK/97bF38f2xXLZ1wjaWmiLQHeCkzsHY0B4YKkG1NERAADhg5F6oV9kI8CVfX2l/kIgv2tJltTKT2hsXGXFjDs/2tvy0JHlNVqJc0wiF1CJTWHIcDH0+nKk5bEB9tfOnSpvA5JEpbiAEBKmB8u2gkYLpaKq2nQWmqkP84WW88RaZKYF2FkWOZj/cL+aF4V0uOkXyj2jrNd8fliaS1SwqQHIgOTQ7DHyta1h3Kq7G7BZ811ncPw6znzO+sAsPNcmaQdVlob2S0CmywsH/rpRBFGdnP8LnD/xGBcKtegpOaPWYbyukYczavCMAf7Chh2nblnYDz+/f1J6PUCStQNWPzjKcx38sJeJpPhxVvTcbpQjQdW7MPyX7Ow9N6+Tl/YA0BKuB8+mDYAKx8YiDfu7iMpKdseP6WnQ7kFRERtxe0ChnfffRfJycnw9vZG//79sXPnTqvnrl+/HmPGjEF4eDgCAgIwZMgQ/PTTT39ib6WRmvQcJHKGQapQETslGXZIkrYW3FhkzRapMwxilddpESrhLrhMZqiSaEtNg07ytma+Sk+bS2YAQ65BgsSlHZEB1pe2GEnZBtcoNcIP52wEDNkVGtHbv7bWPdr2DMPRvGr0jguS3G56XCAOWwkY6rSGINeR9dX9EoOx/1KlxccyThZhtMTEZKNgXy/UNOhM7qwDhu16g1QKh4NyP6UnukT640D2H30WBAGfZ+bhrgHxDrVp9PexXfD8NyfQrBegF4AXvj+FeWO6Ol3j5JZeMUiLDsD0FfuwYP0x/GtyT4S64HeBl6ccT93UDSt+v7CXOhtGRERuFjCsW7cOc+fOxbPPPotDhw7h+uuvx7hx45CTk2Px/B07dmDMmDHYuHEjDhw4gBEjRmDChAk4dOjQn9xzcdSSZxi8RC3HkSpYJSJgcODiM8zfC+V2Ep/LarUI95fWrkrpiTqt7SRaR7Zkk8lgM+E3v0pawrNYl8rqJCd52rtYMy4xk3pRF+qnREWd9aTy04U16CqyCJpZuzYqledKqG/RWqCPAjUNTRbHbW9WOQY5mOSp8JAj1NfLLCgTBAGni2qQFi39NTC6oWuEWTLx2n25uOca53YymXVdMpZuOdey5e63Rwow5PfEWmekRQdgXHo0Zn1yEO+fkmNE13CHZoMsuWdgAj6dNQjLZlzjsoRUIiJynlsFDG+88QZmzZqFBx98EGlpaViyZAni4+Px3nvvWTx/yZIleOqpp3DNNdcgNTUVL730ElJTU/Hdd9/9yT0XR93QBD9v8X/MA1UKVNuZYbC3w40lIb4KmxdzAFBRp5UeMIjYWtUwwyBtKj7SRuKoM8L9bfc3r6LeoYtaX6Unamwsd8ou14gurtaayssDtVYCpzIJNTMuZ2unpGP51ZJqJbTmp7Scf1Lf2Ayv3xONHdEz1vKuRjvOluE6CYWkLjf69+0/WztZqEZqhL9Td9fv6B+HzzPzWoK6mgYdjuRV4RqJSdSXiwzwxs3p0fjH18fxzeF8bDiUj1nXJTvVptHE3jF4597eeKCrHpP7iK/nQERE7slt6jA0NjbiwIEDmD9/vsnxsWPHYteuXaLa0Ov1qKmpQUiI9buMWq0WWu0fF4lqteHCQ6fTQadz/d381pr1zdA3N0Evrko3/LxkuFjSYLNf6nodfL08JPU90NsDZep6m88pqa5HoLenpHaDfTxRVKVBWqT1u+dVmkb4eAqS2g3zU6Cgsg7xQYblC8bnGv9bq22Ct0Iuefwi/ZXILqtpqbR5uUtlNUgJ95XcbmyQElkl1u9KN+iaIBf00On0Fh+3JjHEB+eLqtEjxjyh9UxhFZJCfRx6DyeF+OBMYZXFJOSzRWokBSc71G7f+EDsu1iKG7qEm4xZZq4aveMCHP68DUkOxrbTRegWaRp0nS+pQUKQl8PtXt8pGA+vPoy7+/9xgbxmbzbuGRDn1O8GLzkwODkIa/Zewp394/CfTafxwNAENDWJq4diy+TeUYj0VyCrTIP/3pUOmdAMnYjdysTwkgvw9kCb/14k17n8dyN1fBwz9+NuYya2nzKhvbbDkaigoACxsbH47bffMHTo0JbjL730Ej7++GOcOXPGbhuvvvoqXn75ZZw6dQoREZYT/xYuXIhFixaZHV+9ejVUKtfsVmHNB6fl+Es38ReJBRrgQJkcExKsP6esAdhRKMdtyeLbPVctQ14dMCLG+lvjl3wZkvwFdJKw2cqRchk0TcCQSOvtSn0NAOBgmQx6ARgQbrndIg2wt0SOSUnS2t1XKoNCDvQNtdzu+ktyXBepR4TESYbdxTKoPIHeFtrVC8CyM3LMlvgaAMD+Uhk8ZEC/MPN2dxbJEKIEegRL/7gfLJOhSQAGXvb6CoJhvB5Kk95XAMivM7x/JyaaPn9jrhzdg/RIcnCVT5Me+OiMHHNa9atAA+wplvY5sOSbS3L0DNGjUwBQqwNWnXf852+tWQC+uChHWQOQFiRgVKxb/FomIiI3p9FoMGXKFFRXVyMgwPpFndvMMBhdPvUvCIKo5QBr1qzBwoUL8c0331gNFgBgwYIFmDdvXsvXarUa8fHxGDFiBEJDpRVmkurrioMYP76f6POL1Q3I3nYR48db3/v7SF411OfKMH5EJ9Htni2uwXdHizB+jPXqqkc3ncEtA+KQLGGtfVROFfZmVWD8cOv7qkt9DQAg7FIFjuRVY/zvyy10Oh0yMjIwZswYKBQK7DxfBs8yDcZLqBoMAKFZFTheoMb4a5MsPv7DmsO4d1IvKCUUWAOAkIsVOFloud28ynr0ES5hvAOVUmPzqrHDyljv+uYEpo7oJKkWhVGPCg1W7so261N2hQb9kYPx4x3byaZZL2DvqkMYP76fyZhtWHsMs+/oY7NqtT37mk+g+6Cklt2m3sg4h0cGhaOvg7sZGQ2ua8QTnx/F1Im98e+Np7Hwznj0lVhYzZoJEP/7rCO4/HNGHR/HzP1wzNyPu42ZcSWNPW4TMISFhcHDwwNFRabJgSUlJYiMtL1Dybp16zBr1ix88cUXGD16tM1zlUollErznTkUCkWbD7xcJpf0PcIC5KhpaLb5nJpGPcICfCS1Gx6oQlV9k83nVGqaEBnkK6ndqCAVKjQlVp/TrBfg6eEh+XWODfbDL6fLzZ5nHLOSGh3iQqT1FQDiQ/3wy5kyG/0F/Hyk7+KSHOGPjNOlFtvNrapCpwh/h95rqZGB+GxvrsXnltQ0Ii7Ez6GL0U4RAcivajBr90RhLfomBDv8uVDAUJyptlGAn5ehjTqdIYld5e3c7jh3D0zE2swCPHdLd9Q3NuNYQQ2eGpfm9MV4ZJACT97UDc9+cxITesdgYCfHcyKuFH/G70ZyLY6Z++GYuR93GTOxfXSbpGcvLy/0798fGRkZJsczMjJMlihdbs2aNZgxYwZWr16Nm2++ua276bDGJj08PaRdzHgrPOxuU1pZ14hgifuDB6u8UGkn6VndoEOAhARtwH7Ss9StT40iA7xRXGM96bmgqh6xQdKTk6MCvVFYZbldZ1byRQf6oMBKu5fK6yTN2rQWqFJA3WC+7r2pWQ+5TObwxbJMJoPCQ25WqXvPxXIMSnFu1m3MZYnEW8+UOrxFaWv9EoJRWF2PY3nVeH3zGUwfmuSyO/d9E4LxwbQBmNCbyb5ERHR1cJuAAQDmzZuHjz76CMuXL8epU6fwxBNPICcnB3PmzAFgWE40bdq0lvPXrFmDadOm4fXXX8fgwYNRVFSEoqIiVFdbLxjVXqQWbROrUqNDiEraRbjCQ96yFaMtUi/A7NUgKKuRXtsBAHy8PNBgo92C6gbEOBAwKD09zPbHN6rU6BAk8XU18pDLrAYcF0sdDxgAQC4zBAitXSqvQ0q4420ChoJdmZfVIcivanAoEGttZLcI/Hi8sOXrjceLMLq78wEDALw4OR2f7rmExFAVxrioTSIioquR2yxJAoC7774b5eXleOGFF1BYWIiePXti48aNSExMBAAUFhaa1GR4//330dTUhEceeQSPPPJIy/Hp06dj5cqVf3b3baqud/wC1JbKukaH2m2PVPjSWq3DW3/a4sgsiz25FRrEhzh+sWys8XB5bYicCg3igx1Pru8a5Y/TRTXo2Wqr06N51SZfO2JYl3BsOJTfUnk4t0KDGAmVs63x91agc4Qfdl8sx7lqGZLDfCUVL7Ql2NcL/7mjt0vaIiIiuppJDhjOnTuHXbt2oaioCDKZDJGRkRg6dChSU60nyLrSww8/jIcfftjiY5cHAdu2bWv7DrlI280wNCLY17XtttXGWsXVDYgKdPwi3FbCqKPLUZQKw1Kcyyvu5lY6d2EfHeiD/Kp6xLeq6CwIAgTBPIiQom98MA7lVpkECPsvVeDRkc59PrtF+eNEQXXLa7zxWCHGpUc71abR46NSMXftIRQXy/DpbeKT84mIiOjPIXpJUnV1NSZNmoSuXbti7ty5WL58OT766CPMnTsX3bp1w+TJk0VnWpM5db0OAQ4EDB5ymdVlMwBQpdEh2IEZBrlcZnVZkrqhyeG7wJ4e1vtbpG5AlIN3rYNUlqteOxvcxAT6oKCq3ux4Vmkdkp1Y5tM9JgAnC00/L8VqrcM/v1GfhCAcynH90iGZTIb+iSHYfbEcer2A7WdLMbSTa3YN8/dW4H/39cXsbnr4OVmFmIiIiFxPdMDw2GOPISsrC7t370ZlZSXOnDmDs2fPorKyErt27UJWVhYee+yxtuzrFa2qvtGhGYYgH8sXykbaJvO74+LaVaDKSuJzWa0WoQ4uHQr1U6K81nK7hdUNiHbwgjkyQIkiC9Wey+saJVekbi06yAeF1ebtni+tRadw80JmYvWIMa9GfKpIjW5REgpbWBDmp0RpjbalMnOJugFhTvz8rd0/OAHvbD2P1zPOYFzPKKe2PSUiIiL3Ifov/rfffosPP/wQgwYNMnts0KBBeP/99/HNN9+4tHNXk2qNY0uSglQKVGlcX00wxNf6TknltY0IdSA5GbC9U1JZrdahpGfg952S1ObtFlY5lvBsFB/sg+xyjdnxmoYmp5aQdYn0x5ki04DheF41uluo0ixVemwgjhcYEvt/OlnssoTfCH9v/OPm7kiN8Mf9gxNd0iYRERF1fJJuEdpaB+4uxYY6qup6xy5AA1UKVNfb3gLVESG+XlZnAsprHdv+FADC/bxQWmM5YNDrBXg4uH7fEDCYzwTkVmoQF+x4wNAl0h9ni2tMjulF7CBlj2FLXL3Jkqmj+dVIdzI5GQAm9I7BF5l5EAQBGSeLMbyr62oFpEUHYHLfWH7eiYiIriKiA4YJEyZg9uzZyMzMNHssMzMTc+bMwcSJE13auauJo0nPQT5eVmcYnFm/H2xjhqGszrHtT4Hfl8zYqMXgqMgAJYotLB3KKnNum9KEEBVyK0xnGAqq6xET5PwOQUmhvrhYVgfAULRO16x3aPnY5dKiA6Bu0OG5b45jcEoIVF7MCyAiIiLHiQ4Yli5dipiYGAwcOBAhISHo1q0b0tLSEBISgkGDBiE6OhpvvfVWW/b1ilZdr0OQi5ckaRqb4ePgxWKorxcq6iy3W+5EDkOYv+UlSQ26Znh5Or4mPibIBwXV5snJ2eV1SHKmroFcBv3vuxcZnS6sQddIf4fbNLqucxh+PVcGADicW+WS2QWjxbelY1zPaMwZxl2HiIiIyDmiryaDgoLw448/4vTp09i9ezeKiooAAFFRURgyZAi6devWZp28GqgbHNslKchHYTEpFzBsqRriYP2BYF8vnCq0vOuV0zkMNeYzFyVqLSIDHL9rH+GvtJjDUF7b6PDyKaNwf8OsSIS/oX9H86owygXViId0CsVjaw5h+tAkbDxWiNv6xTrdppHKyxPXdg5zWXtERER09ZJ8+7lbt24MDtqAo+v3A1UKVFtZOlRZ53gxuBCV9RmG0hotwh0OGLwszjAUqR3fIQmwnkMj2HhMrJ6xgTiaW43R3Q39O1VUg0dGdnaqTcBQ+TohRIXPM3NxsbQW3aOdT3gmIiIicjXui+jmglReqLKyrWqlxvEKxyF+1nMYdM16h5cP+Sk9UattMjteUFWPaCdrBQRethVsnbYJKi/ncwIGJYdib1Y5gD9yDZSezrcLAE/e2BV5FRq8MKknE4mJiIioQ2LA4OYM9RJsBAwOLsfxV3qipsH127VauyjOLtcgIcTxyskAkBLmiwuldS1fnyupRecIx2slGKVG+OFMcS0Aw3IkV+Ya+Co9MW9sV5OKz0REREQdCQMGN6fy8oCm0fyOPQBUOFG0TCaTwdImS816wek74d4KOeobm02O5VRokOjkRXOnCD9cKKlt+fpkYQ16xDh/cS+XyxAT6I2LpbXYeKwQY7tHOd0mERERkbtgwNABaJuaHa6aa+vivbRGi3B/x3INrCmv1SLc37kk4vhgFfIqTbcqrdQ0IsjB5VNG6bGBOJZf3fL1qaIalxRCA4AZ1ybhH18fR1aZBulxrpthICIiIuroHNpzs6qqCl9++SUuXLiAJ598EiEhITh48CAiIyMRG+u6nV6uFtX1Oqcvli1xpnIyYLiz3tSsh2erYKZI3dCyW5Cj4kNUyK3UIPWyrUmdnbmIC/ZBbqtApKCqHjFOJFK31i0qAP+a3NOpnZyIiIiI3JHkgOHo0aMYPXo0AgMDcenSJcyePRshISHYsGEDsrOz8cknn7RFP69oageLttlTUadDsIO7JAGGYmglNVrEtEpGLlZrEeXkRXhCiAoXS/9YOtSga4a3wvnJLplMBh+FB+q0TdA0GRKsXZlI3Cnc+XwIIiIiIncj+Spt3rx5mDFjBs6dOwdv7z8uHMeNG4cdO3a4tHNXi0qNDoFOzDDIZDI0680TDgTBsa1ajWKCfFBQZVoMrUjdgMgA55Y5JYf9UeEYAM6X1LrsYvz61HBsO1uGU1UyXNsp1CVtEhEREV3NJAcM+/fvx0MPPWR2PDY2tqWYG0njbHGxAG8F1Fa2VnVGbJAP8i8LGPIr6xEb5FxycmyQD/Iq/2j3dFENukY5XzkZAManR+HzA3nYXSzDuJ5MTiYiIiJyluSAwdvbG2q1eQXgM2fOIDw83CWdutoY6iU4HjAYKxG3JggCnF2NY5hhMK0inVuhQXyIc/US5L/Peuh/nxU5XahGtyjXJCcHqbzwxKjOmJCgh7+3Qyk6RERERNSK5IBh0qRJeOGFF6DTGe5oy2Qy5OTkYP78+bj99ttd3sGrQUVdI0L9HA8YogKUKKo2vbCvrtchwMm8CEtLkup1zVB5OX8hHt8qQfliWR2SQl1Xh6BPfBASXTNhQURERHTVkxwwvPbaaygtLUVERATq6+sxfPhwdO7cGf7+/njxxRfboo9XvPLaRoT4Op4XEBngjWK1acDgii1VI/2VKLqsXVelEA9ICsHeixVo0DXDUy4z2YmJiIiIiDoOybeKAwIC8Ouvv2LLli04ePAg9Ho9+vXrh9GjR7dF/64KlZpGhDixJCkiwBsXLpSZHCut0SLciS1VAcDTQ26STF1dr4Ofi5b5XJ8ahmc2HEOwrxcGpTA5mYiIiKijcvjqb+TIkRg5cqQr+3LVUtfrEODj+IV4VKA3itWmOQwF1Q0m26E6ystD/vu2px64WFqL5DBfp9sEDLkGfkpPvL3lHFY8MNAlbRIRERGR64m6Sn3rrbdEN/j444873JmrlQDnipaF+ylRUmO6dKigqh7DujifhJ4c7oussjqkRQfgdFEN0qJdk5wMAC/emo7qeh1CnNghioiIiIjalqiA4c033zT5urS0FBqNBkFBQQAMlZ9VKhUiIiIYMLQDL085mppN6zAUVNUjJsj5qsSpEX44V1JrCBgK1Zh1XYrTbRopPOROVaImIiIiorYnKtM0Kyur5d+LL76IPn364NSpU6ioqEBFRQVOnTqFfv364V//+ldb9/eK5IpE4svLtpXVahHmRCK1Ua+4IBzKqQQA5FRoEBfs/DInIiIiInIfkremee6557B06VJ07dq15VjXrl3x5ptv4h//+IdLO3c10DQ2wdvLw+l2jLkGRoLwR70DZ3QK98X5klrUNOjg4+XhkjaJiIiIyH1IDhgKCwtbajC01tzcjOLiYpd06mpSotYiwsntTwEgMVSFnApDXQO9/vL5BsfJZDKkhPli0XcnMaZ7pMvaJSIiIiL3IDlgGDVqFGbPno3MzEwIguHCNDMzEw899BC3VnVAsboBkQHO5xokhRmSkwGgoLreJTskGT0+KhWdI/wwoVeMy9okIiIiIvcgOWBYvnw5YmNjMXDgQHh7e0OpVGLQoEGIjo7GRx991BZ9vKIV12gR5YqAIdQXl34PGM6X1KJTuGu2PwWAUD8l5gzvxOJqRERERFchyZv/h4eHY+PGjTh79ixOnz4NQRCQlpaGLl26tEX/rnjF1Q3oEev8VqUp4b74+lA+AEPA0C3KddufEhEREdHVy+FqYV26dGGQ4ALF6gaMTItwup0IfyWKf6/FcLJAjcl9Y51uk4iIiIhIcsAwc+ZMm48vX77c4c5cjYprtC7JYZDJZPBVekLdoENZXSPrGxARERGRS0gOGCorK02+1ul0OH78OKqqqjBy5EiXdexqUadtgp/S4YkeEyO7RuDF708hNcLPJe0REREREUm+Ut2wYYPZMb1ej4cffhgpKa6rAkzS3dwrGqeL1Jh9PceBiIiIiFzDJdveyOVyPPHEE3jzzTdd0dxVo07bBB8XFG0z8lZ44NmbuyPCBUuciIiIiIgAFwUMAHDhwgU0NTW5qrmrQna5BkmhqvbuBhERERGRVZKXJM2bN8/ka0EQUFhYiB9++AHTp093Wceseffdd/Hqq6+isLAQPXr0wJIlS3D99ddbPX/79u2YN28eTpw4gZiYGDz11FOYM2dOm/dTjKyyOiSGuq5eAhERERGRq0kOGA4dOmTytVwuR3h4OF5//XW7Oyg5a926dZg7dy7effddXHvttXj//fcxbtw4nDx5EgkJCWbnZ2VlYfz48Zg9ezY+++wz/Pbbb3j44YcRHh6O22+/vU37KsbxgmpWTyYiIiKiDk1ywLB169a26Icob7zxBmbNmoUHH3wQALBkyRL89NNPeO+997B48WKz8//3v/8hISEBS5YsAQCkpaUhMzMTr732WocIGM4U1eDvY7ijERERERF1XJIDhpEjR2L9+vUICgoyOa5WqzF58mRs2bLFVX0z0djYiAMHDmD+/Pkmx8eOHYtdu3ZZfM7u3bsxduxYk2M33ngjli1bBp1OB4VCYfYcrVYLrVbb8rVarQYAZJXUoFLnCQECBAHQC4blWAIACID+9/8XBLSc0/r/9a3OFQBoGpvh7SmDoG+GTt/s+AtDZnQ6ncl/qePjmLkfjpn74Zi5H46Z+3G3MRPbT8kBw7Zt29DY2Gh2vKGhATt37pTanGhlZWVobm5GZGSkyfHIyEgUFRVZfE5RUZHF85uamlBWVobo6Giz5yxevBiLFi0yO/76N3vh5aOCDGj5B5mV/wcgkxn/XzAcl5k+BgD9QgRs3Jhv92cnx2RkZLR3F0gijpn74Zi5H46Z++GYuR93GTONRiPqPNEBw9GjR1v+/+TJkyYX6c3Nzdi0aRNiY2MldNExMuPV9u8EQTA7Zu98S8eNFixYYJLYrVarER8fjzdnDEdoaKij3aY/kU6nQ0ZGBsaMGWNxFok6Ho6Z++GYuR+OmfvhmLkfdxsz40oae0QHDH369IFMJoNMJrNY0dnHxwdLly4V30OJwsLC4OHhYTabUFJSYjaLYBQVFWXxfE9PT6sX/0qlEkql0uy4QqFwi4GnP3DM3A/HzP1wzNwPx8z9cMzcj7uMmdg+ig4YsrKyIAgCUlJSsG/fPoSHh7c85uXlhYiICHh4uK4I2eW8vLzQv39/ZGRk4NZbb205npGRgUmTJll8zpAhQ/Ddd9+ZHNu8eTMGDBjgFoNIRERERNTeRAcMiYmJAAC9Xt9mnbFn3rx5mDp1KgYMGIAhQ4bggw8+QE5OTktdhQULFiA/Px+ffPIJAGDOnDl4++23MW/ePMyePRu7d+/GsmXLsGbNmnb7GYiIiIiI3ImogOHbb7/FuHHjoFAo8O2339o8d+LEiS7pmCV33303ysvL8cILL6CwsBA9e/bExo0bW4KZwsJC5OTktJyfnJyMjRs34oknnsA777yDmJgYvPXWWx1iS1UiIiIiIncgKmCYPHkyioqKEBERgcmTJ1s9TyaTobm5bbcIffjhh/Hwww9bfGzlypVmx4YPH46DBw+2aZ+IiIiIiK5UogKG1suQ2nNJEhERERER/bnk7d0BIiIiIiLquETNMLz11luiG3z88ccd7gwREREREXUsogKGN998U1RjMpmMAQMRERER0RVEVMCQlZXV1v0gIiIiIqIOyKkcBkEQIAiCq/pCREREREQdjEMBw7Jly9CzZ094e3vD29sbPXv2xEcffeTqvhERERERUTsTXenZ6LnnnsObb76Jxx57DEOGDAEA7N69G0888QQuXbqEf//73y7vJBERERERtQ/JAcN7772HDz/8EPfee2/LsYkTJ6JXr1547LHHGDAQEREREV1BJC9Jam5uxoABA8yO9+/fH01NTS7pFBERERERdQySA4b7778f7733ntnxDz74APfdd59LOkVERERERB2D5CVJgCHpefPmzRg8eDAAYM+ePcjNzcW0adMwb968lvPeeOMN1/SSiIiIiIjaheSA4fjx4+jXrx8A4MKFCwCA8PBwhIeH4/jx4y3nyWQyF3WRiIiIiIjai+SAYevWrW3RDyIiIiIi6oCcKtxGRERERERXNskzDA0NDVi6dCm2bt2KkpIS6PV6k8cPHjzoss4REREREVH7khwwzJw5ExkZGbjjjjswcOBA5ioQEREREV3BJAcMP/zwAzZu3Ihrr722LfpDREREREQdiOQchtjYWPj7+7dFX4iIiIiIqIORHDC8/vrrePrpp5Gdnd0W/SEiIiIiog5E8pKkAQMGoKGhASkpKVCpVFAoFCaPV1RUuKxzRERERETUviQHDPfeey/y8/Px0ksvITIykknPRERERERXMMkBw65du7B792707t27LfpDREREREQdiOQchm7duqG+vr4t+kJERERERB2M5IDh5Zdfxt///nds27YN5eXlUKvVJv+IiIiIiOjKIXlJ0k033QQAGDVqlMlxQRAgk8nQ3Nzsmp4REREREVG7kxwwbN26tS36QUREREREHZDkgGH48OFWHzt8+LAzfSEiIiIiog5Gcg7D5aqrq/Huu++iX79+6N+/vyv6REREREREHYTDAcOWLVtw//33Izo6GkuXLsX48eORmZnpyr4REREREVE7k7QkKS8vDytXrsTy5ctRV1eHu+66CzqdDl999RW6d+/eVn0kIiIiIqJ2InqGYfz48ejevTtOnjyJpUuXoqCgAEuXLm3LvhERERERUTsTPcOwefNmPP744/jrX/+K1NTUtuwTERERERF1EKJnGHbu3ImamhoMGDAAgwYNwttvv43S0tK27BsREREREbUz0QHDkCFD8OGHH6KwsBAPPfQQ1q5di9jYWOj1emRkZKCmpqYt+0lERERERO1A8i5JKpUKM2fOxK+//opjx47h73//O15++WVERERg4sSJbdFHAEBlZSWmTp2KwMBABAYGYurUqaiqqrJ6vk6nw9NPP4309HT4+voiJiYG06ZNQ0FBQZv1kYiIiIjoSuNUHYauXbviP//5D/Ly8rBmzRpX9cmiKVOm4PDhw9i0aRM2bdqEw4cPY+rUqVbP12g0OHjwIJ577jkcPHgQ69evx9mzZ9s0qCEiIiIiutJIrvRsiYeHByZPnozJkye7ojkzp06dwqZNm7Bnzx4MGjQIAPDhhx9iyJAhOHPmDLp27Wr2nMDAQGRkZJgcW7p0KQYOHIicnBwkJCS0SV+JiIiIiK4kLgkY2tru3bsRGBjYEiwAwODBgxEYGIhdu3ZZDBgsqa6uhkwmQ1BQkNVztFottFpty9dqtRqAYYmTTqdz7AegP5VxnDhe7oNj5n44Zu6HY+Z+OGbux93GTGw/3SJgKCoqQkREhNnxiIgIFBUViWqjoaEB8+fPx5QpUxAQEGD1vMWLF2PRokVmx7du3QqVSiW+09TuLp9hoo6PY+Z+OGbuh2Pmfjhm7sddxkyj0Yg6r10DhoULF1q8OG9t//79AACZTGb2mCAIFo9fTqfT4Z577oFer8e7775r89wFCxZg3rx5LV+r1WrEx8djxIgRCA0Ntfu9qP3pdDpkZGRgzJgxUCgU7d0dEoFj5n44Zu6HY+Z+OGbux93GzLiSxp52DRgeffRR3HPPPTbPSUpKwtGjR1FcXGz2WGlpKSIjI20+X6fT4a677kJWVha2bNlic3YBAJRKJZRKpdlxhULhFgNPf+CYuR+OmfvhmLkfjpn74Zi5H3cZM7F9bNeAISwsDGFhYXbPGzJkCKqrq7Fv3z4MHDgQALB3715UV1dj6NChVp9nDBbOnTuHrVu3coaAiIiIiEgip7ZV/bOkpaXhpptuwuzZs7Fnzx7s2bMHs2fPxi233GKS8NytWzds2LABANDU1IQ77rgDmZmZWLVqFZqbm1FUVISioiI0Nja2149CRERERORW3CJgAIBVq1YhPT0dY8eOxdixY9GrVy98+umnJuecOXMG1dXVAIC8vDx8++23yMvLQ58+fRAdHd3yb9euXe3xIxARERERuR232CUJAEJCQvDZZ5/ZPEcQhJb/T0pKMvmaiIiIiIikc5sZBiIiIiIi+vMxYCAiIiIiIqsYMBARERERkVUMGIiIiIiIyCoGDEREREREZBUDBiIiIiIisooBAxERERERWcWAgYiIiIiIrGLAQEREREREVjFgICIiIiIiqxgwEBERERGRVQwYiIiIiIjIKgYMRERERERkFQMGIiIiIiKyigEDERERERFZxYCBiIiIiIisYsBARERERERWMWAgIiIiIiKrGDAQEREREZFVDBiIiIiIiMgqBgxERERERGQVAwYiIiIiIrKKAQMREREREVnFgIGIiIiIiKxiwEBERERERFYxYCAiIiIiIqsYMBARERERkVUMGIiIiIiIyCoGDEREREREZBUDBiIiIiIisooBAxERERERWcWAgYiIiIiIrGLAQEREREREVjFgICIiIiIiqxgwEBERERGRVW4TMFRWVmLq1KkIDAxEYGAgpk6diqqqKtHPf+ihhyCTybBkyZI26yMRERER0ZXGbQKGKVOm4PDhw9i0aRM2bdqEw4cPY+rUqaKe+/XXX2Pv3r2IiYlp414SEREREV1ZPNu7A2KcOnUKmzZtwp49ezBo0CAAwIcffoghQ4bgzJkz6Nq1q9Xn5ufn49FHH8VPP/2Em2+++c/qMhERERHRFcEtZhh2796NwMDAlmABAAYPHozAwEDs2rXL6vP0ej2mTp2KJ598Ej169PgzukpEREREdEVxixmGoqIiREREmB2PiIhAUVGR1ee98sor8PT0xOOPPy76e2m1Wmi12pav1Wo1AECn00Gn00noNbUX4zhxvNwHx8z9cMzcD8fM/XDM3I+7jZnYfrZrwLBw4UIsWrTI5jn79+8HAMhkMrPHBEGweBwADhw4gP/+9784ePCg1XMsWbx4scU+bd26FSqVSnQ71P4yMjLauwskEcfM/XDM3A/HzP1wzNyPu4yZRqMRdZ5MEAShjftiVVlZGcrKymyek5SUhNWrV2PevHlmuyIFBQXhzTffxAMPPGD2vCVLlmDevHmQy/9YddXc3Ay5XI74+HhcunTJ4vezNMMQHx+PwsJChIaGiv/hqN3odDpkZGRgzJgxUCgU7d0dEoFj5n44Zu6HY+Z+OGbux93GTK1WIywsDNXV1QgICLB6XrvOMISFhSEsLMzueUOGDEF1dTX27duHgQMHAgD27t2L6upqDB061OJzpk6ditGjR5scu/HGGzF16lSLAYaRUqmEUqk0O65QKNxi4OkPHDP3wzFzPxwz98Mxcz8cM/fjLmMmto9ukcOQlpaGm266CbNnz8b7778PAPjLX/6CW265xWSHpG7dumHx4sW49dZbERoaajYjoFAoEBUVZXNXJSIiIiIi+oNb7JIEAKtWrUJ6ejrGjh2LsWPHolevXvj0009Nzjlz5gyqq6vbqYdERERERFcet5hhAICQkBB89tlnNs+xl45hLW+BiIiIiIgsc5sZBiIiIiIi+vMxYCAiIiIiIqsYMBARERERkVUMGIiIiIiIyCoGDEREREREZBUDBiIiIiIisooBAxERERERWcWAgYiIiIiIrGLAQEREREREVjFgICIiIiIiqxgwEBERERGRVQwYiIiIiIjIKgYMRERERERkFQMGIiIiIiKyigEDERERERFZxYCBiIiIiIis8mzvDnR0giAAAGpqaqBQKNq5NySGTqeDRqOBWq3mmLkJjpn74Zi5H46Z++GYuR93GzO1Wg3gj+tdaxgw2FFeXg4ASE5ObueeEBERERG5Xk1NDQIDA60+zoDBjpCQEABATk6OzReSOg61Wo34+Hjk5uYiICCgvbtDInDM3A/HzP1wzNwPx8z9uNuYCYKAmpoaxMTE2DyPAYMdcrkhzSMwMNAtBp7+EBAQwDFzMxwz98Mxcz8cM/fDMXM/7jRmYm6IM+mZiIiIiIisYsBARERERERWMWCwQ6lU4vnnn4dSqWzvrpBIHDP3wzFzPxwz98Mxcz8cM/dzpY6ZTLC3jxIREREREV21OMNARERERERWMWAgIiIiIiKrGDAQEREREZFVDBiIiIiIiMgqBgw2vPvuu0hOToa3tzf69++PnTt3tneXyIaFCxdCJpOZ/IuKimrvblErO3bswIQJExATEwOZTIavv/7a5HFBELBw4ULExMTAx8cHN9xwA06cONE+nSUA9sdsxowZZp+7wYMHt09nCYsXL8Y111wDf39/REREYPLkyThz5ozJOfycdSxixoyfs47lvffeQ69evVqKsw0ZMgQ//vhjy+NX4meMAYMV69atw9y5c/Hss8/i0KFDuP766zFu3Djk5OS0d9fIhh49eqCwsLDl37Fjx9q7S9RKXV0devfujbffftvi4//5z3/wxhtv4O2338b+/fsRFRWFMWPGoKam5k/uKRnZGzMAuOmmm0w+dxs3bvwTe0itbd++HY888gj27NmDjIwMNDU1YezYsairq2s5h5+zjkXMmAH8nHUkcXFxePnll5GZmYnMzEyMHDkSkyZNagkKrsjPmEAWDRw4UJgzZ47JsW7dugnz589vpx6RPc8//7zQu3fv9u4GiQRA2LBhQ8vXer1eiIqKEl5++eWWYw0NDUJgYKDwv//9rx16SJe7fMwEQRCmT58uTJo0qV36Q/aVlJQIAITt27cLgsDPmTu4fMwEgZ8zdxAcHCx89NFHV+xnjDMMFjQ2NuLAgQMYO3asyfGxY8di165d7dQrEuPcuXOIiYlBcnIy7rnnHly8eLG9u0QiZWVloaioyORzp1QqMXz4cH7uOrht27YhIiICXbp0wezZs1FSUtLeXaLfVVdXAwBCQkIA8HPmDi4fMyN+zjqm5uZmrF27FnV1dRgyZMgV+xljwGBBWVkZmpubERkZaXI8MjISRUVF7dQrsmfQoEH45JNP8NNPP+HDDz9EUVERhg4divLy8vbuGolg/Gzxc+dexo0bh1WrVmHLli14/fXXsX//fowcORJarba9u3bVEwQB8+bNw3XXXYeePXsC4Oeso7M0ZgA/Zx3RsWPH4OfnB6VSiTlz5mDDhg3o3r37FfsZ82zvDnRkMpnM5GtBEMyOUccxbty4lv9PT0/HkCFD0KlTJ3z88ceYN29eO/aMpODnzr3cfffdLf/fs2dPDBgwAImJifjhhx9w2223tWPP6NFHH8XRo0fx66+/mj3Gz1nHZG3M+DnreLp27YrDhw+jqqoKX331FaZPn47t27e3PH6lfcY4w2BBWFgYPDw8zCLBkpISs4iROi5fX1+kp6fj3Llz7d0VEsG4oxU/d+4tOjoaiYmJ/Ny1s8ceewzffvsttm7diri4uJbj/Jx1XNbGzBJ+ztqfl5cXOnfujAEDBmDx4sXo3bs3/vvf/16xnzEGDBZ4eXmhf//+yMjIMDmekZGBoUOHtlOvSCqtVotTp04hOjq6vbtCIiQnJyMqKsrkc9fY2Ijt27fzc+dGysvLkZuby89dOxEEAY8++ijWr1+PLVu2IDk52eRxfs46HntjZgk/Zx2PIAjQarVX7GeMS5KsmDdvHqZOnYoBAwZgyJAh+OCDD5CTk4M5c+a0d9fIiv/7v//DhAkTkJCQgJKSEvz73/+GWq3G9OnT27tr9Lva2lqcP3++5eusrCwcPnwYISEhSEhIwNy5c/HSSy8hNTUVqampeOmll6BSqTBlypR27PXVzdaYhYSEYOHChbj99tsRHR2NS5cu4ZlnnkFYWBhuvfXWduz11euRRx7B6tWr8c0338Df37/lLmdgYCB8fHwgk8n4Oetg7I1ZbW0tP2cdzDPPPINx48YhPj4eNTU1WLt2LbZt24ZNmzZduZ+xdtufyQ288847QmJiouDl5SX069fPZIsz6njuvvtuITo6WlAoFEJMTIxw2223CSdOnGjvblErW7duFQCY/Zs+fbogCIYtH59//nkhKipKUCqVwrBhw4Rjx461b6evcrbGTKPRCGPHjhXCw8MFhUIhJCQkCNOnTxdycnLau9tXLUtjBUBYsWJFyzn8nHUs9saMn7OOZ+bMmS3Xh+Hh4cKoUaOEzZs3tzx+JX7GZIIgCH9mgEJERERERO6DOQxERERERGQVAwYiIiIiIrKKAQMREREREVnFgIGIiIiIiKxiwEBERERERFYxYCAiIiIiIqsYMBARERERkVUMGIiIiIiIyCoGDERE5LSFCxeiT58+7fb9n3vuOfzlL38Rde7//d//4fHHH2/jHhERXTlY6ZmIiGySyWQ2H58+fTrefvttaLVahIaG/km9+kNxcTFSU1Nx9OhRJCUl2T2/pKQEnTp1wtGjR5GcnNz2HSQicnMMGIiIyKaioqKW/1+3bh3++c9/4syZMy3HfHx8EBgY2B5dAwC89NJL2L59O3766SfRz7n99tvRuXNnvPLKK23YMyKiKwOXJBERkU1RUVEt/wIDAyGTycyOXb4kacaMGZg8eTJeeuklREZGIigoCIsWLUJTUxOefPJJhISEIC4uDsuXLzf5Xvn5+bj77rsRHByM0NBQTJo0CZcuXbLZv7Vr12LixIkmx7788kukp6fDx8cHoaGhGD16NOrq6loenzhxItasWeP0a0NEdDVgwEBERG1iy5YtKCgowI4dO/DGG29g4cKFuOWWWxAcHIy9e/dizpw5mDNnDnJzcwEAGo0GI0aMgJ+fH3bs2IFff/0Vfn5+uOmmm9DY2Gjxe1RWVuL48eMYMGBAy7HCwkLce++9mDlzJk6dOoVt27bhtttuQ+sJ9YEDByI3NxfZ2dlt+yIQEV0BGDAQEVGbCAkJwVtvvYWuXbti5syZ6Nq1KzQaDZ555hmkpqZiwYIF8PLywm+//QbAMFMgl8vx0UcfIT09HWlpaVixYgVycnKwbds2i98jOzsbgiAgJiam5VhhYSGamppw2223ISkpCenp6Xj44Yfh5+fXck5sbCwA2J29ICIiwLO9O0BERFemHj16QC7/475UZGQkevbs2fK1h4cHQkNDUVJSAgA4cOAAzp8/D39/f5N2GhoacOHCBYvfo76+HgDg7e3dcqx3794YNWoU0tPTceONN2Ls2LG44447EBwc3HKOj48PAMOsBhER2caAgYiI2oRCoTD5WiaTWTym1+sBAHq9Hv3798eqVavM2goPD7f4PcLCwgAYliYZz/Hw8EBGRgZ27dqFzZs3Y+nSpXj22Wexd+/ell2RKioqbLZLRER/4JIkIiLqEPr164dz584hIiICnTt3NvlnbRemTp06ISAgACdPnjQ5LpPJcO2112LRokU4dOgQvLy8sGHDhpbHjx8/DoVCgR49erTpz0REdCVgwEBERB3Cfffdh7CwMEyaNAk7d+5EVlYWtm/fjr/97W/Iy8uz+By5XI7Ro0fj119/bTm2d+9evPTSS8jMzEROTg7Wr1+P0tJSpKWltZyzc+dOXH/99S1Lk4iIyDoGDERE1CGoVCrs2LEDCQkJuO2225CWloaZM2eivr4eAQEBVp/3l7/8BWvXrm1Z2hQQEIAdO3Zg/Pjx6NKlC/7xj3/g9ddfx7hx41qes2bNGsyePbvNfyYioisBC7cREZFbEwQBgwcPxty5c3HvvffaPf+HH37Ak08+iaNHj8LTk6l8RET2cIaBiIjcmkwmwwcffICmpiZR59fV1WHFihUMFoiIROIMAxERERERWcUZBiIiIiIisooBAxERERERWcWAgYiIiIiIrGLAQEREREREVjFgICIiIiIiqxgwEBERERGRVQwYiIiIiIjIKgYMRERERERkFQMGIiIiIiKy6v8Bh0vrm0PrNDcAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 900x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m  = 1.0            # system mass in kg\n",
    "fn = 1.0            # natural frequency in Hz\n",
    "wn = 2*np.pi*fn     # natural frequency in rad/s\n",
    "zt = 0.02           # damping ratio to critical (nondim)\n",
    "k  = wn*wn*m        # system stiffness (N/m)\n",
    "\n",
    "print(k)\n",
    "\n",
    "u  = F.sdof_Duhamel(fn, zt, U0=0., V0=0.)/m\n",
    "\n",
    "u.plot_time(fig=1, figsize=(9,3), axis_t=[0, Td, -0.5, 0.5]);\n",
    "\n",
    "print('Maximum response amplitude: {0:6.0f}mm'.format(1000*u.max()))\n",
    "\n",
    "print(1/(2*np.pi*1.))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that ``MRPy`` has calculated the peak response amplitude as 15.4mm.\n",
    "Now we will perform the same calculation by assuming that the applied force function \n",
    "can be approximated by an impulse function at time origin.\n",
    "\n",
    "The impulse function equivalent to the load defined above is given by:\n",
    "\n",
    "$$ F(t) = F_0 t_{\\rm d} \\; \\delta(t) = I \\delta(t) $$\n",
    "\n",
    "where $F_0 t_{\\rm d}$ is the total impulse (integral of the Dirac's delta), which\n",
    "is not unitary in the example. As previously stated, the equivalent initial \n",
    "velocity, $v_0$, is approximatelly:\n",
    "\n",
    "$$ v_0 = \\frac{I}{m} = \\frac{F_0 t_{\\rm d}}{m} = 1{\\rm m/s} $$\n",
    "\n",
    "where $I$ is the given impulse. The response amplitude to this initial condition is:\n",
    "\n",
    "$$ u_{\\rm max} \\approx \\frac{v_0}{\\omega_{\\rm D}} \\approx \\frac{v_0}{\\omega_{\\rm n}} \n",
    "               \\approx \\frac{1.0}{2 \\pi \\cdot 1.0} \\approx 159{\\rm mm} $$\n",
    "\n",
    "which is quite close to the numerical result (the discretization $N$ was intentionally\n",
    "chosen to give the required accuracy).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2. Dynamic factor approach\n",
    "\n",
    "As a rule of thumb, the approximation of a short transient load as an impulse function\n",
    "at time origin provides a good approximation up to $t_{\\rm d} \\leq T_{\\rm n}/4$ (where\n",
    "$ T_{\\rm n}$ is the system natural period of vibration). \n",
    "For instance, by setting $t_{\\rm d} = 0.25$s in the previous example, the numerical\n",
    "result would be $u_{\\rm max} = 35.3$mm, while the impulse function approximation would give\n",
    "$u_{\\rm max} = 39.8$mm, overestimating the correct result (try it by yourself!).\n",
    "\n",
    "To overcome this approximation error, numerically computed results may be provided as \n",
    "design curves, as the one reproduced below (Clough & Penzien, 1995):\n",
    "\n",
    "<img src=\"images/dynamic_factors.jpg\" alt=\"Dynamic factors\" width=\"560px\"/>  <br>\n",
    "\n",
    "The curves provide a _dynamic amplification factor_, to be multiplied by the system \n",
    "static response to the force function maximum amplitude: \n",
    "\n",
    "$$ u_{\\rm max} = A \\frac{F_0}{k} = A \\frac{F_0}{\\omega_{\\rm n}^2 m} $$\n",
    "\n",
    "In the previous example, with $t_{\\rm d} = 0.25$s, the dynamic amplification taken\n",
    "from the respective curve (rectangular function shape) is $A \\approx 1.4$. \n",
    "Hence the peak response amplitude is:\n",
    "\n",
    "$$ u_{\\rm max} = 1.4 \\frac{4.0}{(2 \\pi \\cdot 1.0)^2 \\cdot 1.0} \\approx 142{\\rm mm}$$\n",
    "\n",
    "which is much closer to the correct solution.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.10132118364233778 0.14184965709927289\n"
     ]
    }
   ],
   "source": [
    "k    = ((2*np.pi*fn)**2)*m\n",
    "uest =  F0/k\n",
    "\n",
    "A    =  1.4\n",
    "udyn =  A*uest\n",
    "\n",
    "print(uest, udyn)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Assignments <a name=\"section_5\"></a> \n",
    "\n",
    "1. Utilizar o registro de vibração livre obtido no trabalho anterior.\n",
    "2. Uma vez identificados frequência e amortecimento, calcular:\n",
    "   (a) Máximo deslocamento para uma carga transiente de curta duração, através de          aproximação por resposta impulsiva.\n",
    "   (b) Mesmo cálculo feito por ``MRPy.sdof_Duhamel()``.\n",
    "3. Relatório com descrição do objeto, gráficos e resultados (nome do arquivo       T3_xxxxxxxx.ipynb).\n",
    "\n",
    "Prazo: 20 de maio de 2020.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}