S-Curve motion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# S-Curve motion (3rd order motion) calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### S-Curve calculation for stepper motors\n",
    "This paper shows a way to calculate the complete S-Curve motion without any iteration.\n",
    "\n",
    "S-Curve motion means that, given the speed, acceleration and hyper-acceleration (called jerk), given a starting and a ending point, the calculation results give a smooth motion: at any time the jerk is maintained under a given maximum.\n",
    "\n",
    "**Also, the calculation gives the minimum optimum travel time.** (see README for thanks)\n",
    "\n",
    "During the travel, a S-Curve motion can take up to 7 phases. Depending of the parameters, the phases may be lowered until 4.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
       "    return false;\n",
       "}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
    "    return false;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import math\n",
    "import numpy as np\n",
    "from IPython.display import Image\n",
    "from IPython.core.display import HTML \n",
    "#size of plot\n",
    "plt.rcParams['figure.figsize']=(6,6)\n",
    "plt.rcParams['figure.dpi']=150"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Intro\n",
    "In the remaining part of the document, **the half** of the complete travel is addressed. The other half is exactly the same movement but with jerk inverted.\n",
    "\n",
    "Assume the system can sustain a maximum jerk of $j_{max}$, the motor can have a maximal acceleration of $a_{max}$ and a maximum velocity of $v_{max}$. We suppose also that $x_{start}<x_{end}$. Let define also $x_{end/2}$ the half of the travel requested. \n",
    "\n",
    "There are several motion phases in this half travel movement:\n",
    "\n",
    "The first phase **PH1** is the acceleration motion, from rest, the acceleration raise from $0$ (at $t=0$) to $a_{max}$ at $t=t_0$, during this phase we have $j(t)=j_{max}$.\n",
    "\n",
    "The second phase **PH2** is the motion at maximum acceleration, during this phase we have $j(t)=0$, the end of this phase is noted $t_1$.\n",
    "\n",
    "The third phase **PH3** is the deceleration motion, during this phase we have $j(t)=-j_{max}$, the end of this phase is noted $t_2$.\n",
    "\n",
    "The fourth phase **PH4** is the one at constant speed $v(t)=v_{max}$ ($a(t)=0$) until we reach $x_{end/2}$.\n",
    "\n",
    "This is summarized in the following figure that represents the motion velocity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"IM1.jpg\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Image(url= \"IM1.jpg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculations\n",
    "The calculations must be aware that in some cases, the maximum velocity can't be reached. In the figure above, $x_{end/2}-x_{start}$ can be so small that **PH4** collapse completly. In this case, a new lower $v_{max}$ is calculated in order to make $x_{t2}$ the mid travel point ($x_{end/2}$).\n",
    "\n",
    "Also it is possible that **PH2** collapse completly if the parameters are such that $a_{max}$ cannot be reached before $v_{max}$ is reached. In this case, a new lower $a_{max}$ is calculated in order to have $v_{max}$ at $x_{t0}$.\n",
    "\n",
    "Finally, it's possible that $a_{max}$ can't be reached. In this case a new $a_{max}$ is calculated in order to complete the half travel without **PH2**.\n",
    "\n",
    "This is described here after."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. $x_{end/2}-x_{start}$ big enough to have complete travel at $v_{max}$\n",
    "Here we have all the phases\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**For PH1** the movement equations are:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "j(t) & = j_{max} \\\\\n",
    "a(t) & = j_{max}.t \\\\\n",
    "v(t) & = \\dfrac{j_{max}.t^2}{2}\\\\\n",
    "x(t) & = \\dfrac{j_{max}.t^3}{6}+x_{start}\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "**For PH2**, the change is when $a=a_{max}$ then at $t=t_0=\\dfrac{a_{max}}{j_{max}}$, we have:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "v_{t0} & = \\dfrac{j_{max}.t_0^2}{2}= \\dfrac{a_{max}.t_0}{2}\\\\\n",
    "x_{t0} & = \\dfrac{j_{max}.t_0^3}{6}+x_{start} = \\dfrac{a_{max}.t_0^2}{6}+x_{start}\n",
    "\\end{aligned}\n",
    "$$\n",
    "we note $t'=t-t_0$, the movement equations are:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "j(t') & = 0 \\\\\n",
    "a(t') & = a_{max} \\\\\n",
    "v(t') & = a_{max}.t'+v_{t0}\\\\\n",
    "x(t') & = \\dfrac{a_{max}.t'^2}{2}+v_{t0}.t'+x_{t0}\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "**For PH3**, let $t=t_1$, the change is when $v(t_1)=v_{max}-v_{t0}$ by symmetry. Hence $v_{t1}=v_{max}-\\dfrac{a_{max}^2}{2.j_{max}}$. At that time, we have also $v_{t1}=v_{t0}+a_{max}.(t_1-t_0)$, then we can calculate $t_1$, we have $t_1=\\dfrac{v_{max}}{a_{max}}$. At that point we also calculate $x(t_1)$, we have finally:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "v_{t1} & = a_{max}.(t_1-t_0)+v_{t0}\\\\\n",
    "x_{t1} & = \\dfrac{a_{max}.(t_1-t_0)^2}{2}+v_{t0}.(t_1-t_0)+x_{t0}\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "we note $t''=t-t_1$, the movement equations are:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "j(t'') & = -j_{max} \\\\\n",
    "a(t'') & = -j_{max}.t''+a_{max} \\\\\n",
    "v(t'') & = \\dfrac{-j_{max}.t''^2}{2}+a_{max}.t''+v_{t1}\\\\\n",
    "x(t'') & = \\dfrac{-j_{max}.t''^3}{6}+\\dfrac{a_{max}.t''^2}{2}+v_{t1}.t''+x_{t1}\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "**For PH4**, assume the change is at $t=t_2$, we have $t_2-t_1=t_0$ by symmetry. Hence $t_2=\\dfrac{a_{max}}{j_{max}}+\\dfrac{v_{max}}{a_{max}}$. At $t=t_2$ we have $v(t_2)=v_{max}$, we can calculate $x(t_2)=x_{t2}$:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "v_{t2} & = v_{max}\\\\\n",
    "x_{t2} & = \\dfrac{-j_{max}.t_0^3}{6}+\\dfrac{a_{max}.t_0^2}{2}+v_{t1}.t_0+x_{t1}\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "we note $t'''=t-t_2$, the movement equations are:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "j(t''') & = 0 \\\\\n",
    "a(t''') & = 0 \\\\\n",
    "v(t''') & = v_{max}\\\\\n",
    "x(t''') & = v_{max}.t'''+x_{t2}\\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. $x_{end/2}-x_{start}$ not big enough reach $v_{max}$ during the (half) travel\n",
    "In this case, given the position $x_{end/2}$, we have to calculate a new $v_{max}$ in order to have $v(t_2)=v_{max}$.\n",
    "\n",
    "Let's develop $x_{t2}$ found in the precedings equations, we have $x_{t2}=x_{end/2}$, hence: \n",
    "$$\n",
    "\\begin{aligned}\n",
    "x_{end/2} & = \\dfrac{-j_{max}.t_0^3}{6}+\\dfrac{a_{max}.t_0^2}{2}+v_{t1}.t_0+x_{t1}\\\\\n",
    "& = \\dfrac{-j_{max}.t_0^3}{6}+\\dfrac{a_{max}.t_0^2}{2}+(a_{max}.(t_1-t_0)+v_{t0}).t_0+\\dfrac{a_{max}.(t_1-t_0)^2}{2}+v_{t0}.(t_1-t_0)+x_{t0}\\\\\n",
    "& = \\dfrac{-j_{max}.t_0^3}{6}+\\dfrac{a_{max}.t_0^2}{2}+a_{max}.t_1.t_0-a_{max}.t_0^2+v_{t0}.t_0+\\dfrac{a_{max}.t_1^2}{2}-a_{max}.t_1.t_0+\\dfrac{a_{max}.t_0^2}{2}+v_{t0}.t_1-v_{t0}.t_0+x_{t0}\\\\\n",
    "& = \\dfrac{-j_{max}.t_0^3}{6}+\\dfrac{a_{max}.t_1^2}{2}+v_{t0}.t_1+x_{t0}\\\\\n",
    "& = \\dfrac{-j_{max}.t_0^3}{6}+\\dfrac{a_{max}.t_1^2}{2}+(\\dfrac{j_{max}.t_0^2}{2}).t_1+\\dfrac{j_{max}.t_0^3}{6}+x_{start}\\\\\n",
    "& = \\dfrac{a_{max}.t_1^2}{2}+\\dfrac{j_{max}.t_0^2.t_1}{2}+x_{start}\\\\\n",
    "& = \\dfrac{a_{max}.(\\dfrac{v_{max}}{a_{max}})^2}{2}+\\dfrac{j_{max}.(\\dfrac{a_{max}}{j_{max}})^2.\\dfrac{v_{max}}{a_{max}}}{2}+x_{start}\\\\\n",
    "& =\\dfrac{v_{max}^2}{2.a_{max}}+\\dfrac{a_{max}.v_{max}}{2.j_{max}}+x_{start}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The equation can be rewriten with $v_{max}$ as variable:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "j_{max}.v_{max}^2+a_{max}^2.v_{max}+2.a_{max}.j_{max}.(x_{start}-x_{end/2})=0\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "It's a second order equation, let's calculate the roots:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "v_{max} & = \\frac{-a_{max}^2 \\pm \\sqrt{ a_{max}^4 + 8.a_{max}.j_{max}^2.(x_{end/2}-x_{start})}}{2.j_{max}}\\\\\n",
    "& = \\frac{-a_{max}^2 \\pm \\sqrt{ a_{max}^4.(1 + \\frac{8.j_{max}^2}{a_{max}^3}.(x_{end/2}-x_{start}))}}{2.j_{max}}\\\\\n",
    "& = \\frac{a_{max}^2}{2.j_{max}}.\\left(\\mp \\sqrt{ 1 + \\frac{8.j_{max}^2}{a_{max}^3}.(x_{end/2}-x_{start})}-1\\right)\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "In this equation, we have $x_{end/2}\\ge x_{start}$ then $a_{max}\\ge 0$, hence the root is always positive and $\\ge1$. From the hypothesis in building the equation, $v_{max}$ must be $\\ge 0$, hence the final result is\n",
    "$$\n",
    "\\begin{aligned}\n",
    "v_{max} = \\frac{a_{max}^2}{2.j_{max}}.\\left(\\sqrt{ 1 + \\frac{8.j_{max}^2}{a_{max}^3}.(x_{end/2}-x_{start})}-1\\right)\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "This equation is valid only if **PH2** is still present, in the case of **PH2** is collapsed, we have to recalculate $x_{t2}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. PH2 collapsed\n",
    "When $t_1=t_0$, it remains only the phase **PH1** and **PH3**, let's calculate the value of $x_{t2}$:\n",
    "\n",
    "**For PH1** the motion equations are the same:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "j(t) & = j_{max} \\\\\n",
    "a(t) & = j_{max}.t \\\\\n",
    "v(t) & = \\dfrac{j_{max}.t^2}{2}\\\\\n",
    "x(t) & = \\dfrac{j_{max}.t^3}{6}+x_{start}\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "**For PH3**, the change is when $a=a_{max}$ then at $t=t_0=t_1=\\dfrac{a_{max}}{j_{max}}$, we have:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "v_{t0} & = \\dfrac{j_{max}.t_0^2}{2}= \\dfrac{a_{max}.t_0}{2}\\\\\n",
    "x_{t0} & = \\dfrac{j_{max}.t_0^3}{6}+x_{start} = \\dfrac{a_{max}.t_0^2}{6}+x_{start}\n",
    "\\end{aligned}\n",
    "$$\n",
    "we note $t'=t-t_0(=t-t_1)$, the motion equations are:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "j(t') & = -j_{max} \\\\\n",
    "a(t') & = -j_{max}.t'+a_{max} \\\\\n",
    "v(t') & = \\dfrac{-j_{max}.t'^2}{2}+a_{max}.t'+v_{t0}\\\\\n",
    "x(t') & = \\dfrac{-j_{max}.t'^3}{6}+\\dfrac{a_{max}.t'^2}{2}+v_{t0}.t'+x_{t0}\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "Let's calculate $x_{t2}$, note that by symmetry $t_2-t_0=t_0$\n",
    "$$\n",
    "\\begin{aligned}\n",
    "x(t_0) & = \\dfrac{-j_{max}.t_0^3}{6}+\\dfrac{a_{max}.t_0^2}{2}+v_{t0}.t_0+x_{t0}\\\\\n",
    "& = \\dfrac{-a_{max}.t_0^2}{6}+\\dfrac{a_{max}.t_0^2}{2}+\\dfrac{a_{max}.t_0}{2}.t_0+\\dfrac{a_{max}.t_0^2}{6}+x_{start}\\\\\n",
    "& = a_{max}.t_0^2+x_{start}\\\\\n",
    "& = a_{max}.(\\dfrac{a_{max}}{j_{max}})^2+x_{start}\\\\\n",
    "& = \\dfrac{a_{max}^3}{j_{max}^2}+x_{start}\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "We have $x_{t_0}=x_{end/2}$, we find $a_{max}$:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "a_{max}& = \\sqrt[3]{j_{max}^2.(x_{end/2}-x_{start})}\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "We can calculate the $v_{max}$ also:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "v(t_0)= v_{max}& = \\dfrac{-j_{max}.t_0^2}{2}+a_{max}.t_0+\\dfrac{a_{max}.t_0}{2}\\\\\n",
    " & =a_{max}.t_0\\\\\n",
    "  & =\\frac{a_{max^2}}{j_{max}}\\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's calculate all this in python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "vMprocess = 0.1372\n",
      "jM = 30.0000\n",
      "aM = 2.0801\n",
      "vM = 0.1442\n",
      "t0 = 0.0693\n",
      "t1 = 0.0693\n",
      "t2 = 0.1387\n",
      "v_t0 = 0.0721\n",
      "x_t0 = 0.0017\n",
      "v_t1 = 0.0721\n",
      "x_t1 = 0.0017\n",
      "v_t2 = 0.1442\n",
      "x_t2 = 0.0100\n",
      "t_xend2 = 0.1387\n",
      "Time for a complete travel: = 0.2773\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'j(t)')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAE2CAYAAAATcIcbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeZxkZX3o/8+3q6qruqu6p/eenZ4VBGQdBlEDMyo6ylVcUAkRwcRLEiW/n0lMXMg1xBvvVS9GTfD+DHEjGDMoxoREIgIyogaRGZSRAWFWZqeX2bqqu2v9/v44T3Wfqu6e6enqrqXn+369qussT516vlXV9a1zzvM8R1QVY4wxZq6pq3QFjDHGmNlgCc4YY8ycZAnOGGPMnGQJzhhjzJxkCc4YY8ycZAnOGGPMnGQJrkaJyM0ioiKycoJ1Qbfu9mls93YR0aJl80XkfhE54rb7oUkee5F7fNvpPm+5Tff1mYV6tLjX7JIJ1m0SkU0VqNNyEfmGiOwSkaSI9IrI4yLyP09jG0tF5E4R2S4iIyISF5EnReQ2EZk3m/U3Ji9Y6QqYmvAJ4CrgZuAQsGeSchcBfwl8EzhSjorNAS14r9l+4KmidR8od2VE5CxgC/Ai8Em897obWAtcB/yPKWzjSuB+oBf4W+AZIAS8Avgg0AH88czX3phCluDMVLwMeFpVvzdTGxSRACCqmpmpbc41qvpsBZ7294AY8FpVHfAtv1dE/uxUDxaRVuA+4Dngdaqa8K3+oYh8DnhlqZUUEQFCqpoqdVtm7rJDlGcIEekUkb8XkRdEZEhE9onIt0Rk0Uke0+MOV64Dfssd1lMR6Zmg7M3A193s9uKybvpTIvJREdkNpICXi0hERD4vIs+4w1iHReTfReQc37Yvc49/ywTP+39FpE9EQr5lt4jI0+7QWL+IfHU6h019h4FfKSLfFpFBEXlJRD7m1m8QkV+KSMIdfru06PEiIn8sIs+LSEpEDrnDds351xfY7Yr/g+81u9mtH3eIUkTOFpHvicgxERkWkZ+LyIaiMre77awSke+71/VFEfmEiJzqf74NGAGOFa9Q1dwUXrb3A53AHxUlt/w2Eqr6kKvnOlfPdUX1z7/uPb5le0TkmyLyuyLyG7zPz9vcYfO/KX4eEXmX28bFvmVXicgj7n1MiMiDInJ+0ePeICL/JSLH3ev2vIh8YgpxmypkCa72BcQ75zZ6AwITlMt/cX0M2AD8GbAK+JmIRCbZ9iHgCmAr8Es3fYVbXuz7wF+76XdOUvZm4Brgw+7+IBAGmtxjrwH+EIgAj4vIfABVfRJ4HniP/wlFpB54N7BRVdNu2aeBLwEPA29xcW4A/tPtNU7H3cCvgbcB/wr8LxH5DPB/gM+4OkSBf3V1yvsU8DfAQ8Cbgc+61+D7LtEcAt7uyv5vxl6z709UCRFZCPwUuBC4FXgXXiL6voi8cYKHfA/4EfBWV++/Am46Ray/wNuDu1dErhSR8CnKF7saOKSqm0/zcVOxHvgTvDg2AJuBbwO/PcF7eyPwjKr+EkBErgEeAeJ4n6Mb8D53PxGRJa7McrxDq7vx3tO34L1/0VmIxZSDqtqtBm94X5R6itvtJ3l8AFjiyr3Nt/x272NRUPanwKbTqNPKCdYpXkJrOMU2AkAjMAj8sW/5bcAwMM+37K1uu2vdfA+QBT5RtM1XuXJvLarPpK9PUTyf8C0L4p1bSgPLfMvf4spe5ebbgCTwjaJtvseVe4uvzgq8f4Ln3+R/3YE7gIz/9XWv1/PAU8XvIfC+ou39GvjhKWIW4MtAzm0jCfwE+FMgMoXPwHPA41P8DK9zz7Fukte9x7dsDzAEzJ/kvX2Db1mne3/+3LdsB/BI0WObgX7gC27+Oret5pn4H7Vb5W+2B1f73gZcVnR7xUQFReQP3aG7ON4X5V636uxyVBT4gaoOT1Cvd4nIEyJyzNUrgbcX4a/XN/H29t7pW3Yj8Lyq/sLNX413VOKfivZon8BLmFdOs97/mZ9Q75zhDuAFVd3tK/Mbd7/E3b8CqHf19tuIF+NV06jHlcDPVXWHrz5Z4J+Bi/KHPn2K9wSfAZae7AnU8wfACuCPgO8CK/GS6y9EpAG8c6hFRw7K8V3yc1U9XFTfnwE78T4LedfjPgeurqvw4in+XAwBjzP2ufgVXmLcKCLXiUjXrEZjZp0luNr3jKpu9t/wWsEVEJE/Av4v3qG7t+O1issnwskOUc60cYc2ReTNwL14v/xvAC7HS9J9/nqp6ovAY7gvMhFpwTukeY9vc/kvpB14X1T+WxPQPs16Hy2aT02yDF+d8+f8CmJ2CXLAt/50tBVvzzmMt+fVWrS8uCVrkim+16q6W1XvVNUbgMV4h1dfjtcIBbzDff7XN3+eah9w1lSeYxomih28HxFvFZH8ocQbgR+p6gE3n/9cfJXxn4v/hvtcuB8Ob8D7XrwHOOzOcU7nx4ipAtaK8sxxPd4hmj/NLxCRZWWuw0TXZroe2KGqN+cXuAYjEyWAe/AaY5yF90VUvIeUb/X3esYnIP/6csgnl/nAtvxCt+fQzvS6URxx2ys2H++1nSjmkqlqVkQ+Bfw5cK5b/Pt4PxryDrr7h4GrReRSVR33Q6vIiLuvL1o+2Q+Rya7tdQ9eV4u3i8gTeD+Q/Oca8+/7x1z9io22xFTVR4FH3bnHV+F1lfi+iPSoav9kgZjqZAnuzNEInCha9r4Zfo6ku284jcc04h2y87uRiRvKfAe4E/gd4I3AT9yeXd5DeOeOlqprqVdBP8f74rweb28n7914/3eb3PzpvGY/Bj7kvmz3wGh3i3cDv1TV4vf3tInIAlWdaE8p36r1EICqPj/JJr6C17DnThEp7iaAiDQCr1TVh/H62gGcD/zQV+ya06mzqu4Ukf/C+9ysxjvE/S++Is/jncM7T1U/PcVtJoEfiUgM+DdgGd75OlNDLMGdOX4AfEREPo7XUu41eCfVZ1K+39YHReRuvENAW/XkfZV+gHd46fPAfwBr8M79TNRM/YSI/BteZ+EFwH8vWr/TtW68U0TOxksII3jnxa4GvuJ+oc86VT0iXp+vj4lIAngArz/hX+M12smfH3sJbw/jehHZivflvFsL+6DlfR6vAcZDIvKXeD9YPoD3pX5aSeEkbhORV+KdK8yfk7oAb+9tgLGuIBNycb8DrzXiUyLyd4x19F4L/AFeP7mHVfWQiPwY7zXqx2u88x5g+TTqfQ9e69mXA99T1bivTioiHwT+zbVy/TZesurG65O3V1X/RkT+AO983AN4h1o78Pb6DroYTK2pdCsXu03vxslbLAYpaiWIt4fw/+Gd2xrESybLJih3O9NsRenK/iVwAK8142hLODf91xOUr8P70j+Id9L/x8DFeL+4vzFB+WvctgpaVBaVuRFvDyqB1yz8Obw9v8W+MqfTinJl0fJNwE+LlvVQ1BoS77zYH+PtQaTw9n6+RFErPbzWoM/iJRMFbvY9z6aismfjNfk/jpe8fw5sKCpzu9tOsGj5N4A9p4j5crxWlM/g/chI4zVG+gaw4jQ+n2e513wn3l5qHHgS+Kg/frzze//unusw8L/w+tJN1Irymyd5vlb3PAq8fpIyV+B97o+6124PXiK/wrf+3/CSW9K9X98Bzq7U/7ndSruJe2ONMcaYOcVaURpjjJmTLMEZY4yZkyzBGVPlROT3ReQLJ1n/ZhG5t5x1MqYW2Dk4Y6qYa/W3E3iFqh6QsQGaQ+q7EoOIPAPcoKpbK1JRY6qQ7cEZU92uBX6jY6NyTOafgVvKUB9jaoYlOGMqTERWuMu+XOLmF4p3CaB1eB3af+wr/pi7P+Yu53KFm9/EzPWFM2ZOsARnTIWp6k7gI8A33UgfXwfuVtVNeB2X/aOG5AcGblHVmKo+7uafA3omGHDZmDOWjWRiTBVQ1X9wA08/gbucjlvVgtcx/1TyZVoYPySbMWck24Mzpnr8A964jH+n3liI4I260TT5Q0bly4wb4syYM5UlOGOqgBvU9wt4l3S5XUTyV1PYijfWZN5kzZ5fhjcMl+29GeNYgjOmOnwR2Kyq78cbiPnLbvkDFF4ctQ/vignFAxJfhe/CrMYYS3DGVJyIXAtsAP7QLfoT4BIR+R28gYjPEZGFAKo6BHwK+JmIHBOR/EVrfxv4+/LW3JjqZh29jalyInILcK6qfmiS9W8GblTVd5W3ZsZUN0twxhhj5iQ7RGmMMWZOmlKCE5ENIvK8iOwQkY9OsD4sIve69U+48fIQkXYRedSNuHDnzFbdGGOMmdwpE5yIBPCuQvxG4Fzgt0Xk3KJivwccVdWVwOeBz7jlI8D/AD48YzU2xhhjpmAqI5msBXao6i4AEdmINwDss74y1wK3u+n7gDtFRFQ1AfxURFZOtUIdHR3a09Mz1eKTSiQSRKPRkrdTLSye6mbxVDeLpzpkckoynSWZyZHNKZ1NYaC0eLZs2dKvqp0TrZtKglsE7PPN7wcun6yMqmZE5DjQDvRPpYKuldgtAN3d3dxxxx1TedhJxeNxYrFYydupFhZPdbN4qpvFUz45VQaGlUOJHAfjysFEjkPxHAcTORLpsXLz6uH/rG9EREqKZ/369S9Otq4qxqJU1buAuwDWrFmj69atK3mbmzZtYia2Uy0snupm8VQ3i2fmpTI5XhxIsLMvzo5e77a9N86uvgTD6exoubZoPSs7m1mzOsbKrrHbguYIdXUCzF48U0lwB4AlvvnFbtlEZfaLSBCYBwzMSA2NMcZUzImRNDt64+zsjbOzL8GO3ji7+uK8eGSIbG6sm9milgZWdMW4fFl7QSJri9ZXrO5TSXBPAqtEZBleIrseuKGozP3ATcDjwHXAj9Q62BljTE1QVQ4dH2Fnn5fIdvTF2dmbYEdfnL7B5Gi5UEBY1hHl7PlNXHPBAlZ0xljRGWN5Z5RouCoOCBY4ZY3cObVbgQeBAPA1Vd0mIp/EGzvvfrwBYu8RkR3AEbwkCICI7AGagXoReSvwelV9tvh5jDHGzK7iw4r+PbJEauywYlMkyMquGOtWd7Kiy0tiK7tiLGltIBione7TU0q5qvoA3qCv/mWf8E2PAO+c5LE9JdTPGGPMaZrqYcWF8yKs6IrxzjVLWNEVY2VnjBVdUTpjYUSkghHMjOrbpzTGGHNKp3NYsafdO6z4ppcvYGVXdR9WnElzOzpjjKlxiWSG3f0JdvUn2NUXZ3d/gqd3DdP3owcnPKx41erO0SS2ojPK0rbGmjqsOJMswRljTIVlc8qBo8Ps7Pea2ecT2a6+BIdPjIyWE4GF8xpoCcG685e482NRVnbF5sxhxZlkCc4YY8rkSCLF7n7vvNiuvgS7XUJ7cWCIVDY3Wq45EmR5Z4xXrmhneWeU5e6QYk97lEgo4PqNnVfBSGqDJThjjJlBI+kse48MsauvKJH1Jzg2NDaURyggLG1rZHlnjNec0zWayJZ1RGmP1tve2AywBGeMMacplcmx7+gQLw4k2N2fv0+wZyDB/qPD+HsBdzWFWd4Z5U0vX8DyjqiXyDpiLK6xJve1yBKcMcZMIJnJsu/IMHtc4npxYIg9A970gaPD+Frb0xQO0tMR5aIlrbz94sWjSWxZZ5TYHG+pWM3slTfGnLFG0ln2Hx0a3QvbM5BgT7+XyA4eK0xizZEgyzqiXLyklbddtIiejihntUdZ1hGltTFkhxSrkCU4Y8ycdnwozZ7jWR749SH2Hhli7xGXzPqHOHi88HDivIYQPR1RLj2rlXdcspiejkYvibVHabEkVnMswRljaloqk+PgseHR5LXvyBD7jnrTeweGODGS8Qo+/hQArY0hzmqPcllPKz0di+lpj9LTEaWnvZGWxsoNDGxmniU4Y0xVU1UGEqmx5OUSmTc/zKHjhYcS6wN1LG5rYGlbI5csbWVJayMnDu3ijb+1liVtDTRFQpULxpSVJThjTEXlckrvYJIDx4bYf3SYA8eGOVB0P+QbsQOgsynM0rZG1i5rY0lbI0vdbUlbA91NY9cZy9u0aS/nLmwuZ1imCliCM8bMqmQmy6FjI6MJa7+7P3jMS16Hjg+TzhZeXaulMcSilgaWdUR59aqO0QS2tK2Rxa2NNNQHKhSNqSWW4Iwx05bJ5ugdTHLo+AiHj49w+MQIh48Pc+j4WELriycLGnKIQHdThEWtDVy0pIVrLljAopYGFrU2sLilgYUtDXN+EGBTHvYpMsZMaCSd5aUTY4krn8QOHR8eXdY3mCw4/wUQCdUxvznCwpYGrlrdyaLWBl8Ca2T+vAj1QevgbGafJThjzjDDqSx9g0l6B0fcfbJg/qUTSfb1Jxj8wQ/GPbYpEmR+c4T58yKcPb+J+fMaWDDPm18wL8L85gjzGqw5vakOluCMmQMy2RxHhlIcSaS8ZHUiSV/cu88nrvxtMJkZ9/g6gY5YmK7mMPPnRegODnHJOctd4mpgvktiNiqHqSX2aTWmCiUzWY4kUgzEvaR1JJFiIJHiSCI54fLjw+kJtxOtD9DVHKEzFuZlC5u50iWxzlh4dHlnU5i2aD0BX8tDb7T6VeUK15hZYQnOmFmiqgynsxwfTnNsKD16f2I4zbFhLyn51+WnjyRSxCfYywII1AmtjfW0R+tpi9bzsoXNo9Pe/VgC62wKW2MNc0azT78xE1BVhlJZ4skMgyMZ4skM8ZEM8WS6aD7DYDLDrr0j/OOeJ12SSnF8OMPx4dS45u9+gTphXkNo9NbaWM+yjmhBsmqL1tMeG0tgzZHQuD5expiJTSnBicgG4ItAAPiKqn66aH0Y+EfgUmAAeLeq7nHrPgb8HpAF/h9VfXDGam/OeNmct5c0nPJuQ+kMQ6ksI6ksQ6ksQ+kswylvWb7ckLuNpLMMuXWJ5FjSGkxmSCQz41oHTiQSqiMWDhHI5eisG2FeQ4hz5jfT7JJWS6O7d/PNvmWxcNAaYxgzi06Z4EQkAHwJuBrYDzwpIver6rO+Yr8HHFXVlSJyPfAZ4N0ici5wPXAesBB4WERWq2rhsASm6uVySjqXI5NVMtmx6XQ2RyanZLI5Ulm3PpcjXVQuvz6ZzpHM5nh2T5ptj+4gmcmRzGRJZXLedLpoPpMdXZ7Kuvm0t24k7a07HSLQGArQUB+ksT5AY32ASChALBykqylCLBIkFg7SNHofIhYJ0hQOjq7Lr4+Gg4Tc9by8c1a/NRsvvTFmmqayB7cW2KGquwBEZCNwLeBPcNcCt7vp+4A7xftpei2wUVWTwG4R2eG29/jMVH9i2w4e5/6dKbZmtwOMdjJV1DftuAVaODthWf86/wrFO6RVWGbix/g7vOpJnjs/n9/2gYNJ/qPvaXI5JadKViGnSi6nZHNKzs1n3fqxaS85ZdU3XVRGFbJuOptTLzkVJbDsVHZnTtdvngegPlhHOFhHOBhw93XeslCAcKCOWDhIe3Rs/Wj5UICGUGA0UeWTVkP9BMtD3vJwsM72mow5Q0wlwS0C9vnm9wOXT1ZGVTMichxod8t/XvTYRcVPICK3ALcAdHd3s2nTpilWf2I/3p/mX7anYfsLUypf/HWX//6TicoUrfMvH7esaHv+TUy0bLLnERFEc9T1H6DOPU+deDcRqAPqRLxpNz86XVDeK1OP77EBf3khIHi3uvx9gIBAsM6/XIqWjT1u3LI6CLr5oEAoAKE6ITWcYF5TlGCd97zj5dxtChRIupsz4m5HpraFksXj8ZI/t9XE4qluFs/UVEUjE1W9C7gLYM2aNbpu3bqStndlTnn1pk1cddVVo7/WRxOIUJO/4L1DYOsqXY0ZY/FUN4unulk8UzOVBHcAWOKbX+yWTVRmv4gEgXl4jU2m8tgCW7Zs6ReRF6dQr1PpAPpnYDvVwuKpbhZPdbN4qlsp8Zw12YqpJLgngVUisgwvOV0P3FBU5n7gJrxza9cBP1JVFZH7gW+JyN/gNTJZBfziZE+mqp1TqNMpichmVV0zE9uqBhZPdbN4qpvFU91mK55TJjh3Tu1W4EG8bgJfU9VtIvJJYLOq3g98FbjHNSI5gpcEceW+jdcgJQN80FpQGmOMKYcpnYNT1QeAB4qWfcI3PQK8c5LHfgr4VAl1NMYYY07bXL5mxV2VrsAMs3iqm8VT3Sye6jYr8YjqLPRvMsYYYypsLu/BGWOMOYNZgjPGGDMn1VyCE5ENIvK8iOwQkY9OsD4sIve69U+ISI9v3cfc8udF5A3lrPdkphuPiFwtIltE5Nfu/jXlrvtkSnmP3PqlIhIXkQ+Xq84nU+Jn7gIReVxEtrn3KlLOuk+khM9cSETudnE85wZSr7gpxHOliDwlIhkRua5o3U0ist3dbipfrSc33XhE5CLfZ22riLy7vDWfWCnvj1vfLCL7ReTO035yVa2ZG143hZ3AcrwRp54Gzi0q8wHgy276euBeN32uKx8GlrntBGo4nouBhW76fOBApd+fUmPyrb8P+A7w4VqOB6+V8lbgQjffXuOfuRvwxpYFaAT2AD01EE8PcAHeFU+u8y1vA3a5+1Y33VrD8awGVrnphcAhoKVW4/Gt/yLwLeDO033+WtuDGx34WVVTQH7gZ79rgbvd9H3Aa0UKB35W1d1AfuDnSpp2PKr6S1U96JZvAxrEu2xRpZXyHiEibwV248VUDUqJ5/XAVlV9GkBVB7Ty/UBLiUeBqHijFTUAKeBEeao9qVPGo6p7VHUr4wc3fQPwkKoeUdWjwEPAhnJU+iSmHY+qvqCq2930QaAXmJGBM0pQyvuDiFwKdAM/nM6T11qCm2jg5+LBmwsGfgb8Az+f6rHlVko8fu8AnlLvqg2VNu2YRCQGfAT4qzLUc6pKeY9WAyoiD7pDMH9ehvqeSinx3Ack8PYM9gJ3qGq5xrOeTCn/17X6nXBKIrIWb49p5wzVa7qmHY+I1AGfA6Z9qqIqBls20yci5+Fdf+/1la7LDLgd+LyqxqUGB8SeQBB4NXAZMAQ8IiJbVPWRylZr2tbiXbh4Id4hvZ+IyMPqLqVlqoOILADuAW5S1dO7YGJ1+QDwgKrun+73Qa3twZ3OwM9IiQM/l0Ep8SAii4HvAe9V1Ur/UssrJabLgc+KyB7gQ8DHxRsmrpJKiWc/8Jiq9qvqEN5oQJfMeo1PrpR4bgB+oKppVe0FfgZUejzEUv6va/U7YVIi0gx8H7hNVX9+qvJlUEo8VwC3uu+DO4D3isinT+vZK3kCchonLIN4J4KXMXbC8ryiMh+k8AT5t930eRQ2MtlF5U/4lxJPiyv/9kq/LzMVU1GZ26mORialvEetwFN4DTKCwMPANTUcz0eAr7vpKN4YsxdUezy+st9gfCOT3e59anXTbTUcTz3wCPChSsYwU/EUrbuZaTQyqfgLMI0X7E3AC3jHlm9zyz4JvMVNR/Ba4O3Au3LBct9jb3OPex54Y6VjKSUe4C/wzof8ynfrqnQ8pb5Hvm3cThUkuBn4zL0Hr8HMM8BnKx1LiZ+5mFu+DS+5/VmlY5liPJfh7U0n8PZEt/ke+7suzh3A+yodSynxuM9auug74aJajadoGzczjQRnQ3UZY4yZk2rtHJwxxhgzJZbgjDHGzEmW4IwxxsxJluCMMcbMSZbgjDHGzEmW4IwxxsxJluCMMcbMSZbgjDHGzEmW4IypciLy+yLyhZOsf7OI3FvOOhlTC2wkE2OqmIjkL3nyClU94K6uvRsIqXcpm3y5Z4Ab1LuuljEG24MzptpdC/xGVU81Avs/A7eUoT7G1AxLcMZUmIh8RETuK1r2RRH5W+CNwI99qx5z98dEJC4iV7j5TcA1s15ZY2qIJThjKm8j8CYRaQIQkQDwLuBbwMvxrn6Rd6W7b1HVmKo+7uafA3rc9cCMMViCM6biVPVFvOvGvc0teg0wpN4FK1uAwSlsJl+mZeZraExtsgRnTHX4FvDbbvoGNw9wFGiawuPzZY7NcL2MqVmW4IypDt8B1onIYrw9uXyC2wqs9pWbrNnzy4A9qnpi9qpoTG2xBGdMFVDVPryGIl8Hdqvqc27VA8BVvqJ9QA5YXrSJq4D/nOVqGlNTLMEZUz2+BbyOsb03gH8HzhGRhQCqOgR8CviZiBwTkVe4cr8N/H05K2tMtbOO3sZUORG5BThXVT80yfo3Azeq6rvKWzNjqpslOGOMMXOSHaI0xhgzJ1mCM8YYMydZgjPGGDMnBStdgWIdHR3a09NT8nYSiQTRaLT0ClUJi6e6WTzVzeKZeapep0zFmyhuzTHWvEPH/hYVqg96+1ilxLNly5Z+Ve2caF3VJbienh42b95c8nY2bdrEunXrSq9QlbB4qpvFU93mSjyqSjKT46FHH+OiNZczlMoylMowlMoynMqSzORIZrKMpL37ZCZHMp1jJJMl6Vs2ks6XzZF005lcjkxWSWdzZHNKOqtkcr7pbI50Tsm6W6k6YvVs/ourgdLeHxF5cbJ1VZfgjDFmLlFVRtI5BkfSnBjJcGIkzeBIxpsf9u7z817CypJIZRhKZhlKe/f5+UQqw2hu+dGjU65DKCBEggHCoTrCwQDhYB31wToiIW+6KRKkPlBHoE4IBeoIBsSbrvOmg3VCMOCbrqsjFBAC7r5OhDqBujpBABFBBARvuYhbBtS5dZFQYDZe7gKW4IwxZgpSmRzHhlIcHUpzdChVNJ3m2FCKE8NFCczdp7Mn3+OpE4iFg8TCQRrDQaL1ARrrg3Q3RWjsGJuPhgM01Ac4+OJuLjzvHKLhIA31ARpD3novgY0lrnAwQH3QS1xnIktwxpgzTi6nHB9OM5BI0jeYYiCRpH8wyZFEUdIaTnE04SWvRCo76fbCwTpaG+tpbgjSHAnREatnWUeU5oYgTZEQTRFveVMkSHNDiOaItzy/rLE+gMjUk9CmTftZt2bJTLwUc5olOGPMnJDO5uiPJxmIp+hz9958kmd2jvDVnU/QH08xEPcSWWaC80giMK8hRGtjPS2NIbqaIqzubqK1sZ7WxhAtjfWj61oaQ255PQ31s3+4zZw+S3DGmKqWzSkDiSS9J5K8dGKEl9x97+DY9Esnkgwkkkw0MFMkVEcsqCwOZFjUEuHCxfNoj9XTEQvTHgvT4aY7YmHmNYTO2MN5c1FJCU5ENgBfBALAV1T100XrrwS+AFwAXK+q95XyfMaYuSWVyXH4+AgHjg1z6PgwB66BKA0AACAASURBVI8Nc9glrF533xdPjmu1JwLt0TDdzWG6myNcsHgeXU0RuprDdLrE5d3XEw0HXSu9V1UoSlMp005wIhIAvgRcDewHnhSR+1X1WV+xvcDNwIdLqaQxpvbkcspAIsXBY17iOnh8ZNx0f3z8XldLY4hul6xWdzfR3RyhuzlMV3NkdLojFiYUsHEqzMmVsge3FtihqrsARGQjcC0wmuBUdY9blyvheYwxVSiXU3oHk+w7OsS+I0PsOzLMvqNDHDg6zMHjwxw6NkIqW/iv3xAKsLAlwsKWBs45u4sFbnrhvAYWtkRYMK/BzmeZGVNKglsE7PPN7wcuL606xphqcnwozd4jQ2NJ7OhYItt/dJhUpjCBdTeHWdzayAWLW9hwfoRFLQ0scMlr4bwGWhpDp9Va0JhSTPtyOSJyHbBBVd/v5m8ELlfVWyco+w3gPyY7B+eud3ULQHd396UbN26cVp384vE4sVis5O1UC4unutVyPIMp5XAi527K4aEcL8UzDCSF4Uxh2WgIOhrq6GwQOhrq6GoUOhqEzoY62huE+kB1Jq9afn8mYvGMWb9+/RZVXTPRulL24A4A/o4Yi92y06aqdwF3AaxZs0ZnYkiduTI0T57FU92qPZ6RdJY9Awl29yXY1Z9gV1+CXf1xdvcnODaUHi0XCghL2hppaxzhNRcsZmlbI4tbG1nS1sCStkaaI6EKRjF91f7+nC6LZ2pKSXBPAqtEZBleYrseuGFGamWMOW25nHLw+DC7XQLb3Z9gZ5+XxA4cGy5ozNHdHGZ5R4xrXr6AZR1RVnTGWNYRZXFrA8FAnfvCOb9ywRgzA6ad4FQ1IyK3Ag/idRP4mqpuE5FPAptV9X4RuQz4HtAKvFlE/kpVz5uRmhtzhjo+lGZXf3w0ieWn9wwkGEmPnROL1gdY3hnj0rNaue7SxSzvjLG8I8qyjijRsHWBNXNfSZ9yVX0AeKBo2Sd800/iHbo0xpyGZCbLviND7OzL742NJbSBRGq0XKBOWNrWyPKOKK9e2cFytye2ojNKZ1PYGnSYM5r9jDOmQnI55fCJEbcXlj8/5h1S3HdkCH/f5o5YmOWdUa4+t5vlnVGWdcRY3hllaVuj9QczZhKW4IyZRapeZ+fdLoHtdg099gyMP6TYEAqwrCPK+Yvmce2FC0f3xpZ1Rmu2cYcxlWQJzpgZkEgrT+875iUyd8u3WhxMjrW1D9YJS9sbWdbuHVLs6YiyvCNKT0eU+c0R6mwcRGNmjCU4Y6Ygnc1x6NgIe48Mjd7yHZ/3Hhnymto/8jPAGydxUUsDyzqivO2SRSxzCWx5R5RFLV4rRWPM7LMEZwzeocTjw+lxCSw/ffDYSMGAv6GAuP5hjVyweB7po4d57dqXs6wjypK2xrJcrdgYc3KW4MwZIZPN0TuY5OCxYQ6428Fjw964ice8gX/9hxIBOmL1LG5t5OIlrVx7YSNL27yEtrS9kfnNkYLLqmzaNMC68+aXOyxjzElYgjNzwuBIevSyK/nkdfDYCAeOevOHT4yMu+RKS2OIRS0NLG1v5IoV7SxubWCpS2BLWhutr5gxNc7+g01VS2dz9A0mOXxihN4TIxw+PsJhd5HLw8dH3MUuR0iksgWPC9QJ85u9wX7XLmtjYUuERS2N7r6BhS0NlsCMmePsP9xURCqr7D86RH88Rf+gd1HL/sEkLw2OcPi4S2AnRia8XlgoIHQ1RZg/L8I5C5q46uxO5jd78/nk1V10CNEYc+axBGdmzHAqy0AiSX88Rd9gkn6XtPrj+QSW8qYHk975roceHbeN1saQu6hlhHMXNNM9L+KSV3g0qbU11ltzemPMKVmCMxMaSWc5kkhxdCjF0USaI0Mpjo7OpzgylC6aTxV0WvZrjgTpaArTGQvzsoXNXBkLc6LvAJdfcA4dMe/qzJ1NYdpj9YSD1vrQGDMzLMHNUbmcEk9lODGc5vhwmhPDGU6MpDkxnObESMbdT7BuOM3RoTTD6eyk257XEKK1MURrtJ75zRHOmd9MW9Sbb2usp7PJS1odTWE6Jklamzb1se6ypbP5EhhjznAlJTgR2QB8Ee9qAl9R1U8XrQ8D/whcCgwA71bVPaU851ykqiQzOZLpHCOZLMOpLIlUhkQyf5/hqf1pdv10N4lkhkQq6+69dUOpLPFkhqGkdz84kmYwmRl37qpYUyRIcyREc0OI5khw9Hpf+WTV2ujd2qL1owmtpSFkHZWNMTVh2glORALAl4Crgf3AkyJyv6o+6yv2e8BRVV0pItcDnwHeXUqFpyKRzNA3lGPvwBA5VXfzEklOIaeKFt0Xl8nkcmSySjqbI51VMrnc6HQ6W7jOm8+RzinpTI5MThlJZ93NS1qj0xMsT2Zyp0xGADzjvbT1gTqi4QCN9UFi4SDRcIBYOEh3U4TGcKAgaXn3IZobvGQ2z83HIkFrhGGMmdNK2YNbC+xQ1V0AIrIRuBbwJ7hrgdvd9H3AnSIiqlP6Op+2/9h6kI88NgyPjW/EMJvqA3UEA0KwToiEAu5W590HAzRFgnQ2hd18XeH6UICwW9YQChANe8mr0SWvrU89yeuuejWN9UHqg7YHZYwxpyLTzTUich2wQVXf7+ZvBC5X1Vt9ZZ5xZfa7+Z2uTH/Rtm4BbgHo7u6+dOPGjdOqU97hRI5nDg/REPGuhyV44wPWAQgIUOfuZdy8jM4H6yBQBwHxBskN5JeJ188qKN76oLjHz+K1t+LxOLFYbNa2X24WT3WzeKqbxTNm/fr1W1R1zUTrqqKRiareBdwFICJ969evf3EGNtsB9J+yVO2weKqbxVPdLJ7qVko8Z022opQEdwBY4ptf7JZNVGa/iASBeXiNTSalqp0l1GmUiGyeLKvXIounulk81c3iqW6zFU8pJ3OeBFaJyDIRqQeuB+4vKnM/cJObvg740WyffzPGGGOghD04Vc2IyK3Ag3jdBL6mqttE5JPAZlW9H/gqcI+I7ACO4CVBY4wxZtaVdA5OVR8AHiha9gnf9AjwzlKeowR3Veh5Z4vFU90snupm8VS3WYln2q0ojTHGmGpmHaqMMcbMSZbgjDHGzEk1l+BEZIOIPC8iO0TkoxOsD4vIvW79EyLS41v3Mbf8eRF5QznrPZnpxiMiV4vIFhH5tbt/TbnrPplS3iO3fqmIxEXkw+Wq88mU+Jm7QEQeF5Ft7r2KlLPuEynhMxcSkbtdHM+JyMfKXfeJTCGeK0XkKRHJuAEq/OtuEpHt7nZT8WMrYbrxiMhFvs/aVhGZ9WERp6KU98etbxaR/SJy52k/uarWzA2vteZOYDlQDzwNnFtU5gPAl9309cC9bvpcVz4MLHPbCdRwPBcDC930+cCBSr8/pcbkW38f8B3gw7UcD14jrq3AhW6+vcY/czcAG910I7AH6KmBeHqAC/AGfr/Ot7wN2OXuW910aw3HsxpY5aYXAoeAllqNx7f+i8C3gDtP9/lrbQ9udPxLVU0B+fEv/a4F7nbT9wGvFW8MrWvx/jmTqrob2OG2V0nTjkdVf6mqB93ybUCDeFdvqLRS3iNE5K3AbryYqkEp8bwe2KqqTwOo6oCqTn4dovIoJR4FouIN2tAApIAT5an2pE4Zj6ruUdWtQPEFC98APKSqR1T1KPAQsKEclT6Jacejqi+o6nY3fRDoBWZk4IwSlPL+ICKXAt3AD6fz5LWW4BYB+3zz+92yCcuoagY4jvfLeSqPLbdS4vF7B/CUqiZnqZ6nY9oxiUgM+AjwV2Wo51SV8h6tBlREHnSHYP68DPU9lVLiuQ9I4O0Z7AXuUNUjs13hUyjl/7pWvxNOSUTW4u0x7Zyhek3XtOMRkTrgc8C0T1VUxViUZvpE5Dy8yxC9vtJ1mQG3A59X1bjM4sDVZRQEXg1cBgwBj4jIFlV9pLLVmra1QBbv8Fcr8BMReVjdFUVMdRCRBcA9wE2qOm6vqIZ8AHhAVfdP9/ug1vbgTmf8S6Rw/MupPLbcSokHEVkMfA94r6pW+pdaXikxXQ58VkT2AB8CPi7eaDmVVEo8+4HHVLVfVYfwBkW4ZNZrfHKlxHMD8ANVTatqL/AzoNLjIZbyf12r3wmTEpFm4PvAbar68xmu23SUEs8VwK3u++AO4L0i8umTP6RIJU9ATuOEZRDvRPAyxk5YnldU5oMUniD/tps+j8JGJruo/An/UuJpceXfXun3ZaZiKipzO9XRyKSU96gVeAqvQUYQeBi4pobj+QjwdTcdxbv24wXVHo+v7DcY38hkt3ufWt10Ww3HUw88AnyokjHMVDxF625mGo1MKv4CTOMFexPwAt6x5dvcsk8Cb3HTEbwWeDuAXwDLfY+9zT3ueeCNlY6llHiAv8A7H/Ir362r0vGU+h75tnE7VZDgZuAz9x68BjPPAJ+tdCwlfuZibvk2vOT2Z5WOZYrxXIa3N53A2xPd5nvs77o4dwDvq3QspcTjPmvpou+Ei2o1nqJt3Mw0EpwN1WWMMWZOqrVzcMYYY8yUWIIzxhgzJ1mCM8YYMydZgjPGGDMnWYIzxhgzJ1mCM8YYMydZgjPGGDMnWYIzxhgzJ1mCM6bKiUiniPxGRBomWR926yt9aRRjqoolOGOq30eBb6jqMICIbBKR9+dXqneZpK+5csYYxxKcMVXMXcT2JuCbpyj6LeCmKrnorTFVwRKcMVVARD4qIjtFZFBEnhWRt7lVlwPHVHW/K/cp4LeAO0UkLiJ3Arj1R4FXVKL+xlQjS3DGVIedeIlrHt4Vzb/pLlz5cryrXwCgqrcBPwFuVdWYqvqvl/cccGH5qmxMdbMEZ0wVUNXvqOpBVc2p6r3AdrwraLcAg1PczKArb4zBEpwxVUFE3isivxKRYyJyDDgf6MA77Ng0xc00Acdmq47G1BpLcMZUmIicBfwDcCvQrqoteBdIFWArsLroIZNdxPFleFdMNsZgCc6YahDFS1p9ACLyPrw9OPCuqN0iIot85V8Clvs34Na3AT+f9doaUyMswRlTYar6LPA54HG85PVy4GduXQr4BvAe30O+CFwnIkdF5G/dshuAu12fOGMMIKqTHe0wxlQDN0LJT4CL8529i9aH8Q5NXqmqveWunzHVyhKcMcaYOckOURpjjJmTLMEZY4yZkyzBGWOMmZNmPcGJyBIRedSNr7dNRP7f2X5OY4wxZtYbmbjx9Bao6lMi0gRsAd7qmkaP09HRoT09PSU/byKRIBqNlrydamHxVLeZjEeBVCZHMp1lJJMjmcmSziihoBAOBogE6wiHAtQH65AZecbx7P2pbhbPmC1btvSr6oTXQgyWVKspUNVDwCE3PSgizwGLgAkTXE9PD5s3by75eTdt2sS6detK3k61sHiq23TiGUln2dkXZ0dvnJ29cXa46T39Q6SyOcA7xLKsOcKi1gYOHB3m8IkRUkAKqA/U0dPRyMquGCs7Y6zoirGyK8aKzhiRUKDs8VQzi6e6lRKPiLw46bpydhMQkR7gMeB8VT3hW34LcAtAd3f3pRs3biz5ueLxOLFYrOTtVAuLp7qdLJ5EWjkUz3EwkeNgXDmUyHEwnqN/WEfH3BKgq1FYEK1jYayOhTFhYbSOBbE6GoJj+2nDGW9bB+I5DiWUg/EchxI5eocKt9XRICyM1bntjW0rGpraPt+Z9P7UIotnzPr167eo6pqJ1pUtwYlIDPgx8ClV/ZfJyq1Zs0ZtD248i6e6Pfroo5x36RXs8O2J7eiNs703Tt/g2OAi9cE6lndEvb0u362nPVrSXtdIOsuegcTo8+Zvu/oTpDK50XKdTWFWdsZY1R0b3fNb2RWjsymMyFjym2vvj8VT3Urcg5s0wc36IUpXgRDwXeCfTpbcjKl22Zxy4OgwO/oG2f5SfDSh/ebgEMMPPjJarikcZEVXjKtWdxYkkiVtjQTqZv7MWSQU4Jz5zZwzv3lcffcdGRqXeL/31AEGk5mx+kaCBfUc7s2wfGCIRa0Ns1JfY8ph1hOceD8Lvwo8p6p/M9vPZ8xMSGay7OkfGtsbcslhV1+cpG+PqCMWZmVXlCsWBLnyotWs7GpiVXeMrqI9okoJ1Ak9HVF6OqK8ju7R5apK72DSJenB0fgefb6P72zZD8AXnnqUcLCO5Z2Fe3sru2L0dDQSDpZ2ns+Y2VaOPbhXATcCvxaRX7llH1fVB8rw3MacVDyZ8Rp4+JLYzt44Lx4ZIpvzDt+LwKKWBlZ2xXj1yvaxQ4udTcxrDAHuEMurllUylNMiInQ3R+hujvDqVR0F644NpfjOgz+hadGq0dfll3uP8u9PHxwtE6gTzmprHG3Ykk9+K7pixMJlOTBkzCmVoxXlT2HWWjMbMyUD8eS4w3Q7euMcOj4yWibo9nbOnt/ENRcsGG2RuKIzRkP9mbO30tJYz6rWAOvWLi1YPpzyWn3uLDrH+Ohvesnkxs7lL5gXGX3tVnWPJb/2WLjcoZgznP3UMnOGqnLw+Ij3xfvSYMEX8dGh9Gi5xvoAKzpjvGJ5++gX8cquGGe1NxIK2OA+k2moD3D+onmcv2hewfJ0NseLA97hXP9r/u3N+xhKZUfLtTaGRvd+veTXxMquGAvnRaricK6ZeyzBmZqTyeZ48cgQ218q/ELd2Ref8At1w/nzR5PYqu4mFjRHqLOGEzMmFKgbTVx+uZxy6MQI218aLEh+P3jm8IQ/OPzJz35wmJlgCc5UrYkOie3ojbNnIEE6O/6Q2LsvW1JwPsgOiVVWXZ2wqKWBRS0NrDu7q2DdRIeMf75rgO/98sBomVBA6Gkv7FJxJh4yNtNnCc5U3PGhNDv6BguS2PbeOAeODZPvplkncJb7snvdud3WqKHGtcfCtMfCXL68vWD5RI1+fnN4kAe3HSZ/ms/f6KewT99Yox9jwBKcKZN8s/T8+bHHtiX58guPs6M3QX98rCN0vln6xUtbeeelS6xZ+hkmFg5y4ZIWLlzSUrB8sm4bj+8cmLDbxtievNdtwy7sfGayBGdmVDan7D/qnR/zH37a2Rsv6FjcEIRzFuZYf3bnaBJb1dVkHYvNhMLBAGfPb+Ls+U0Fy/Md77f3DhYkv3/71UEGRwo/b2c/+zNfF4/Z7XhvqoMlODMtyUyW3f1THxrqrRcvKmgyvm3L46xf/6oKRmDmgkCdsLS9kaXtjbz2ZYUd2fsGx87zbfrl8wwHA/z4hT7ucx3ZYfKh05Z1RO2IwRxgCc6c1OBImp19xYlskL1HhgrOiSxubWBVVxNXru4cG9m+MzbpOZFnrVm4mUUiQldzhK7mCK9c2cHS5B7WrXsFAMeH0+Ou4LB1/3G+/+tDBed8l7a5KzV0NfkauURpith5vlphCc6gqgwkUuP2xnb0xjl8YqwjdCggLOuIcu7CZt5y4cIZvTyLMeUyryHEpWe1culZrQXLJ7t80Y9f6CtotTu/OTKW8HwNXdqj9dafr8pYgjuD5HLKwePDbM//A/vOWRzz9UuK1gdY0RXjlSvaR5PYqq4YS9saCVq/JDNHRUIBzls4j/MWFnZkz2Rz7HUDVm/3Jb/vbN5HwtfvsqUxVDBeZz75LWppsH6XFWIJbg7yRpYoOqzYF2dnb4Lh9Ng/ZFu0npWdMd54/gJW+c4/LLCRJYwZFQx4LXuXd8Z4/Xljy1WVQ27kHH/y++GzL7HxyX2j5RpCAZZ3Rgv+x7yO7FHryD7LLMHVsKFUhl1F58e29w7y4sBQwdiAC+dFWNndxNq17QX/YG3R+grW3pjaJiIsbGlgYUsDV67uLFh3pPiQf1+cJ/cc5V9/NTZgdbBOOKu9cbQFcf7/cnlnlMZ6+2qeCfYq1oBjQym2H81y+Bd72e77pzlwbHi0TCD/z9IZ4w3nzS8Y+SFqHaGNKau2aD1rl7WxdllbwfJEMjNuZJ7tvXEefq539OoV4OvI7r91xmi1H6Wnxb75qoSq8tKJZGF/Hjd+X3885Ur9mkiojuUdMS49q3VsaCp3Rej6oB3uMKaaRcNBLljcwgWLCzuypzJjpxX8P2Kf2D3ASNrfkb2eFZ0xGjNJdod2j/7/z2+20woTsQRXZtmcjp6wLjw/Fifu6wjd7K6w/JpzuljV1cTQS7t4+2tfaSesjZmD6oN1rOpuYlV3E2/0Lc/llAPHhsd9XzxxKMOj+54dLRdzV5D3N3JZ6RqGnckd2S3BzZKRdFFHaJfEdvUlSGXHfpF1NYVZ1R3jHZcsGmt51RWjM1Z4RehNm/aypK2xEqEYYyqkrk5Y0tbIkrZG1p8zNmD1o48+ynlrrhjr0uC+Y366o4/vPuXryB6oY9kkHdnPhK49luBKNDiSLhwfzw1Rta+oI/TSNu/82FWrOwv6j81rsE6jxpjTIyJ0NUXoaorwyhWFV2Q/4ftOyie/Zw4e54FnCjuyL3HfSf4f1iu7YjTPoY7sluCmQFXpj6fY3js4bqTzl06MDRSc/7V0/sJ5XHvRotETw8s7z4xfS8aYymuOhLhkaSuXLB3fkT1/VMnfF/Yn2/sLjip1N4cLxuuc7KhSLbAE5zPZ8e4dvXGODxd2hF7ZFeNVKzsKmvguaW2wjtDGmKoUCQV42YJmXraguWB5Jptj39Hx33vfferAhO0C/F0aVnZVd0f2siQ4EdkAfBEIAF9R1U+X43kn42+x5E9iO/viBS2W2qP1rOiK8d8uWFDwhlqLJWPMXBF0R56WdUS5+tzCAasPnxgpTHy9cR75zUvcu3msI3u+ZffYVUHGOrJXumX3rCc4EQkAXwKuBvYDT4rI/ar67MkfWbqhVIadvYnRi2nmL+Gyt6gjdL7PySuW+zpCW58TY8wZTERYMK+BBfMa+K1VhR3ZjyZSBZfD2tEbZ8uLR7n/6bGO7P6+ucVXZS9X39xyPMtaYIeq7gIQkY3AtcCsJbj/2HqQv9w0xMAPHhxdlh81YFVXjDeeP3/0CsArumzUAGOMOR2t0Xoui7ZxWU9hR3b/6Er+Pr0/+k3vuJ2Kly1o4h/eu2ZWj4bJbF/pVkSuAzao6vvd/I3A5ap6q6/MLcAtAN3d3Zdu3LixpOd8/kiWh3YPs7SlngXROhbF6uhqFIJVepx4KuLxOLFYrNLVmDEWT3WzeKpbrcWTySm9Q8rBeI6DiRyH4jnSObj14ghQWjzr16/foqprJlpXFbsuqnoXcBeAiPStX7/+xRnYbAfQPwPbqRYWT3WzeKqbxVOFvjs2WUo8Z022ohwJ7gCwxDe/2C2bkKp2TrbudIjI5smyei2yeKqbxVPdLJ7qNlvxlKOJy5PAKhFZJiL1wPXA/WV4XmOMMWewWd+DU9WMiNwKPIjXTeBrqrpttp/XGGPMma0s5+BU9QHggXI8l89dZX6+2WbxVDeLp7pZPNVtVuKZ9VaUxhhjTCXYuFLGGGPmJEtwxhhj5qSaS3AiskFEnheRHSLy0QnWh0XkXrf+CRHp8a37mFv+vIi8oZz1nsx04xGRq0Vki4j82t2/ptx1n0wp75Fbv1RE4iLy4XLV+WRK/MxdICKPi8g2915Fyln3iZTwmQuJyN0ujudE5GPlrvtEphDPlSLylIhk3MAT/nU3ich2d7upfLWe3HTjEZGLfJ+1rSLy7vLWfGKlvD9ufbOI7BeRO0/7yVW1Zm54rTB3AsuBeuBp4NyiMh8AvuymrwfuddPnuvJhYJnbTqCG47kYWOimzwcOVPr9KTUm3/r7gO8AH67lePAacW0FLnTz7TX+mbsB2OimG4E9QE8NxNMDXAD8I3Cdb3kbsMvdt7rp1hqOZzWwyk0vBA4BLbUaj2/9F4FvAXee7vPX2h7c6LiWqpoC8uNa+l0L3O2m7wNeK95gZ9fi/XMmVXU3sMNtr5KmHY+q/lJV8yObbgMaRCRcllqfXCnvESLyVmA3XkzVoJR4Xg9sVdWnAVR1QFWzZar3ZEqJR4GoiASBBiAFnChPtSd1ynhUdY+qbgVyRY99A/CQqh5R1aPAQ8CGclT6JKYdj6q+oKrb3fRBoBeYkYEzSlDK+4OIXAp0Az+czpPXWoJbBOzzze93yyYso6oZ4DjeL+epPLbcSonH7x3AU6qapPKmHZOIxICPAH9VhnpOVSnv0WpAReRBdwjmz8tQ31MpJZ77gATensFe4A5VPTLbFT6FUv6va/U74ZREZC3eHtPOGarXdE07HhGpAz4HTPtURVWMRWmmT0TOAz6Dt7dQ624HPq+qcZkb19sLAq8GLgOGgEdEZIuqPlLZak3bWiCLd/irFfiJiDys7kohpjqIyALgHuAmVR23V1RDPgA8oKr7p/t9UGt7cFMZ13K0jDuUMg8YmOJjy62UeBCRxcD3gPeqaqV/qeWVEtPlwGdFZA/wIeDj4o2CU0mlxLMfeExV+1V1CG+wg0tmvcYnV0o8NwA/UNW0qvYCPwMqPR5iKf/XtfqdMCkRaQa+D9ymqj+f4bpNRynxXAHc6r4P7gDeKyKnd7HsSp6AnMYJyyDeieBljJ2wPK+ozAcpPEH+bTd9HoWNTHZR+RP+pcTT4sq/vdLvy0zFVFTmdqqjkUkp71Er8BReg4wg8DBwTQ3H8xHg6246indNxwuqPR5f2W8wvpHJbvc+tbrpthqOpx54BPhQJWOYqXiK1t3MNBqZVPwFmMYL9ibgBbxjy7e5ZZ8E3uKmI3gt8HYAvwCW+x57m3vc88AbKx1LKfEAf4F3PuRXvltXpeMp9T3ybeN2qiDBzcBn7j14DWaeAT5b6VhK/MzF3PJteMntzyodyxTjuQxvbzqBtye6zffY33Vx7gDeV+lYSonHfdbSRd8JF9VqPEXbuJlpJDgbqssYY8ycVGvn4IwxxpgpsQRnjDFmTrIEZ4wxZk6yBGeMMWZOsgRnjDFmTrIEZ4wxZk6yBGeMMWZOsgRnjDFmTrIEZ4wxZk6yBGdMFXFXY17nps8Vkc0yyVDqItLtrqxdDdcBNKbqWIIzpoqo6nmqusnN/k+8a64pgIjsEZHX+cq+BDwK3FL2ihpTAyzBGVOF3DW91gP/eoqi/wT8Mk4jfAAADApJREFU/uzXyJjaYwnOmCri20u7Gu8q7SNu+T3AUuDfRSTuuzr4E8ByETmrMjU2pnpZgjOmOr0c77JOAKjqjcBe4M2qGlPVz7rlGbxLvVxYkVoaU8UswRlTnVqAwSmWHXTljTE+luCMqU5HgaYplm0Cjs1iXYypSZbgjKlOW4HVRcvGXZ1YRILASuDpclTKmFpiCc6Y6vQQcImIRHzLXgKWF5VbC+xR1RfLVjNjaoQlOGOqkOvj9iPgWt/i/w38hYgcE5EPu2W/A3y53PUzphaI60NqjKkCIrIXeI+qPiYi5wJ3A2t1gn9UEekCfgxcnO9OYIwZYwnOmCohIp14XQHOVtW9la6PMbXODlEaUwVE5DJgO/B3ltyMmRm2B2eMMWZOsj04Y4wxc1Kw0hUo1tHRoT09PSVtI5nOMTCYIByeO1cRSSaTFk8Vm+14ovVBIqHy/R5NJBJEo9GyPd9ss3iqWynxbNmypV9VOydaNyMJTkSWAP8IdON1Rr1LVb8oIrcD/x3oc0U/rqoPnGxbPT09bN68uaT63PvkXj7y3V+P7xVbw+qZoJdvDbN4Ts+lqzu5+3fXzuIzFNq0aRPr1q0r2/PNNounupUSj4hM2gd0pvbgMsCfqupTItIEbBGRh9y6z6vqHTP0PFNy7UWLaDi6k1e98pXlfNpZ9bP/+i+Lp4rNZjy33LOFkXR2VrZtzFw2IwlOVQ8Bh9z0oIg8ByyaiW1PRyQUoLleaI/NnUNgFk91m814GkIBhi3BGXPaZrwVpYj0AI8B5wN/AtwMnAA24+3lHZ3gMbfgrkrc3d196caNG0uuRzweJxaLlbydamHxVLfZjOdzm0dIpJVPXNEwK9ufiL0/1c3iGbN+/fotqrpmwpWqOmM3IAZsAd7u5rv///buP0aOso7j+Ptzu3fXaxHaCiI/bRtbTUtqtaVIJUBEaEEjKiQgEVp/JCL2D/5AfpUokhgVsYYEE1LCLxVTtIhIbATaEEERkLZQetS211KxQCwUA4G21929r3/MU7ocd727nZ3ZmbnvK9l0dmZ2nufT2d3nnpnZZ4AS0dWaPwLuGGobs2fPtmZ49NFHm7KdrPA82ZZknm/e9bSdc/NjiW1/IL5/ss3zHAA8Y4O0J027LEtSO3AfcI+Z/SE0nv81s5qZ9QG3EQ0M65wbgXJbG9VakS7JcS4dTWngJAm4HdhoZkvr5h9Vt9qXgQ3NKM+50aRUEtW+vlZXw7ncadZVlJ8BLgael/RsmHct8FVJs4iuoN4OfLtJ5Tk3arS3iWqf9+CcG6lmXUX5N0ADLDrob96cc0Mr+SFK5xriQ3U5l3HtfojSuYZ4A+dcxpXa5D045xrgDZxzGddeavNzcM41IJUGTtICSZsk9Ui6Oo0ynSuKcpuo1vwQpXMjlXgDJ6kE/BI4G5hOdGXl9KTLda4oop8JeA/OuZFKowc3F+gxs21mtg9YDpybQrnOFUJ7mx+idMXy45UbOXPpXxMvJ/E7eks6H1hgZt8Kzy8GTjKzxXXr+FiUQ/A82ZZknvu37OOBrRXunD+WaEyF5Pn+yba857n9+V66d9VYevpYILmxKDNxw1MzWwYsA5gzZ4414z5Hfr+kbPM8w7ehbwts3cwpp55Geymd68J8/2Rb3vPc9+o6Dut9890MSeVJ49PyMnBc3fNjwzzn3DCUQ6PmPxVwRbFnX43OcvLNTxoN3D+BqZImS+oALgT+lEK5zhVCuS06LOk/9nZF0Vut0dVRSrycxA9RmllV0mLgIaJb59xhZt1Jl+tcUbzbwHkPzhXE3kqNMeUCNHAAZrYSH5fSuYaUwiHKivfgXEHsrfRx+CHJNz8+kolzGdceenA1/6mAK4i9lRpj2pPvwXkD51zG+UUmrmj2VGp0eQPnnDtwkYk3cK4Y9lb66MxDAyfpZ5L+JWm9pPsljQ/zJ0naI+nZ8Lg1fnWdG33Kpf0Xmfg5OFcMvZUaY9rz8TOBR4ATzGwmsBm4pm7ZVjObFR6XNqEs50ad/T24ih+idAWxt5qTc3Bm9rCZVcPTJ4l+yO2ca5JyW/Qx9YtMXBFUa31UapbKObimjkUp6UHgXjP7jaRJQDdRr+4t4Doze3yQ1/lYlEPwPNmWZJ71r1VZuqaX6z49ho+OT/5LAXz/ZF2e8+ypGt9ZtZsLPtbB2ZPbgRaPRSlpFfDhARYtMbMHwjpLgCpwT1j2KnC8me2SNBv4o6QZZvZW/434WJRD8zzZlmSe8pbXYc1T7Ow4mq6usYmU0d/mlzYz7YjJqZSVBs+THbt7a8BGZnx8KqefPAlI7vMzrAbOzD53sOWSFgFfAM6w0CU0s16gN0yvkbQVmAY8E6fCzo02Rx7aSZtg2WPb0i34hQ3plpc0z5Mpx4zvSryM2D8ll7QAuBI4zcx2180/AnjDzGqSpgBTgZQ/oc7l39QjP8C6759Fb6WWWplPPPEE8+bNS628pHmebCmX2pg4riPxcmKfg5PUA3QCu8KsJ83sUknnATcAFaAP+IGZPTiM7b0G/DtWpSKHA683YTtZ4XmyzfNkm+fJtjh5PmJmRwy0IPEbnraKpGcGO/GYR54n2zxPtnmebEsqj49k4pxzrpC8gXPOOVdIRW7glrW6Ak3mebLN82Sb58m2RPIU9hycc8650a3IPTjnnHOjmDdwzjnnCil3DZykBZI2SeqRdPUAyzsl3RuWPxXGxNy/7Jowf5Ok+WnWezCN5pF0pqQ1kp4P/3427boPJs4+CsuPl/S2pCvSqvPBxHzPzZT0D0ndYV+NSbPuA4nxnmuXdHfIsVHSNf1f2wrDyHOqpLWSqpLO77dsoaQt4bEwvVoPrtE8kmbVvdfWS7og3ZoPLM7+CcsPlbRD0i0jLtzMcvMASsBWYArQATwHTO+3zmXArWH6QqLBnwGmh/U7gclhO6Uc5/kkcHSYPgF4udX7J26muuUrgN8DV+Q5D9FIQeuBT4TnH8z5e+4iYHmYHgtsByblIM8kYCbwK+D8uvkTiUZXmghMCNMTcpxnGjA1TB9NNB7w+LzmqVt+M/Bb4JaRlp+3HtxcoMfMtpnZPmA5cG6/dc4F7g7TK4AzJCnMX25mvWb2ItATttdKDecxs3Vm9kqY3w10SepMpdYHF2cfIelLwItEmbIgTp6zgPVm9hyAme0ys/TG2xpYnDwGjJNUBrqAfUR3CmmlIfOY2XYzW080olK9+cAjZvaGmf2P6N6WC9Ko9EE0nMfMNpvZljD9CrATGHCEjxTF2T8oGqj/SODhRgrPWwN3DPCfuuc7wrwB17HoPnVvEv3lPJzXpi1OnnrnAWstGuC61RrOJOkQ4CrghynUc7ji7KNpgEl6KByCuTKF+g4lTp4VwDtEPYOXgJvM7I2kKzyEOJ/rvH4nDEnSXKIe09Ym1atRDeeR1Ab8HGj4VEXswZZda0maAfyUqLeQd9cDvzCzt0OHLu/KwCnAicBuYLWkNWa2urXVathcoEZ0+GsC8LikVWbmg6hniKSjgF8DC83sfb2iHLkMWGlmOxr9PshbD+5l4Li658eGeQOuEw6lHEY0EPRwXpu2OHmQdCxwP3CJmbX6L7X94mQ6CbhR0nbgcuBaSYuTrvAQ4uTZATxmZq9bdKeNlcCnEq/xwcXJcxHwFzOrmNlO4O9Aq8dDjPO5zut3wqAkHQr8mehenU82uW6NiJPnZGBx+D64CbhE0k9GVHorT0A2cMKyTHQieDIHTljO6LfOd3nvCfLfhekZvPcik220/oR/nDzjw/pfafV+aVamfutcTzYuMomzjyYAa4kuyCgDq4DP5zjPVcCdYXoc8AIwM+t56ta9i/dfZPJi2E8TwvTEHOfpAFYDl7cyQ7Py9Fu2iAYuMmn5f0AD/2HnAJuJji0vCfNuAL4YpscQXYHXAzwNTKl77ZLwuk3A2a3OEicPcB3R+ZBn6x4fanWeuPuobhvXk4EGrgnvua8RXTCzAbix1VlivucOCfO7iRq377U6yzDznEjUm36HqCfaXffab4ScPcDXW50lTp7wXqv0+06Yldc8/baxiAYaOB+qyznnXCHl7Rycc845NyzewDnnnCskb+Ccc84VkjdwzjnnCskbOOecc4XkDZxzzrlC8gbOOedcIf0f16kCnKBzQwUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "########################################################\n",
    "# Edit this:\n",
    "########################################################\n",
    "############## system/motor parameters (absolute values)  ################\n",
    "jM = 30     # m/s^3\n",
    "aM = 3     # m/s^2\n",
    "vM = 1  # m/s\n",
    "\n",
    "startX = 0\n",
    "endX2 = 0.01\n",
    "\n",
    "########################################################\n",
    "############### Calculations\n",
    "########################################################\n",
    "# 1. One time calculations\n",
    "t0 = aM / jM\n",
    "t1 = vM / aM\n",
    "\n",
    "t2 = t0 + t1\n",
    "\n",
    "vt0 = jM*t0*t0/2\n",
    "xt0 = jM*t0*t0*t0/6 + startX\n",
    "\n",
    "vt1 = aM*(t1-t0) + vt0\n",
    "xt1 = aM*(t1-t0)*(t1-t0)/2 + vt0*(t1-t0) + xt0\n",
    "\n",
    "vt2 = vM\n",
    "xt2 = -jM*t0*t0*t0/6 + aM*t0*t0/2 + vt1*t0 + xt1\n",
    "\n",
    "#################\n",
    "# Let's detect if half travel (endX2) is >= xt2\n",
    "if endX2 >= xt2:\n",
    "    # Yes, let's calculate endtime via x(t) = vM*(t-t2) + xt2\n",
    "    endtime = t2 + (endX2-xt2) / vM\n",
    "else:\n",
    "    # No, let's calculate new vM:\n",
    "    vM = aM*aM/2/jM*(math.sqrt(1+8*jM*jM/(aM*aM*aM)*(endX2-startX))-1)\n",
    "    print(\"vMprocess = {:.4f}\".format( vM ))\n",
    "    # and recalculate the affected variables\n",
    "    t0 = aM / jM\n",
    "    t1 = vM / aM\n",
    "    # If t1<t0, then we have to lower aM in order to have at least t1=t0\n",
    "    if t1 < t0:\n",
    "        aM =  ( jM*jM*(endX2-startX) )**(1. / 3.0)\n",
    "        vM = aM*aM/jM\n",
    "        #recalculate t0/t1\n",
    "        t0 = aM / jM\n",
    "        t1 = t0    \n",
    "\n",
    "    t2 = t0 + t1\n",
    "    \n",
    "    vt0 = jM*t0*t0/2\n",
    "    xt0 = jM*t0*t0*t0/6 + startX\n",
    "\n",
    "    vt1 = aM*(t1-t0) + vt0\n",
    "    xt1 = aM*(t1-t0)*(t1-t0)/2 + vt0*(t1-t0) + xt0\n",
    "\n",
    "    vt2 = vM\n",
    "    xt2 = -jM*t0*t0*t0/6 + aM*t0*t0/2 + vt1*t0 + xt1\n",
    "\n",
    "    endtime = t2 + (endX2-xt2) / vM\n",
    "    \n",
    "print(\"jM = {:.4f}\".format( jM )) \n",
    "print(\"aM = {:.4f}\".format( aM ))    \n",
    "print(\"vM = {:.4f}\".format( vM ))\n",
    "print(\"t0 = {:.4f}\".format( t0 )) \n",
    "print(\"t1 = {:.4f}\".format( t1 ))    \n",
    "print(\"t2 = {:.4f}\".format( t2 ))\n",
    "\n",
    "print(\"v_t0 = {:.4f}\".format( vt0 )) \n",
    "print(\"x_t0 = {:.4f}\".format( xt0 ))\n",
    "print(\"v_t1 = {:.4f}\".format( vt1 )) \n",
    "print(\"x_t1 = {:.4f}\".format( xt1 ))\n",
    "print(\"v_t2 = {:.4f}\".format( vt2 )) \n",
    "print(\"x_t2 = {:.4f}\".format( xt2 ))\n",
    "print(\"t_xend2 = {:.4f}\".format( endtime ))\n",
    "print(\"Time for a complete travel: = {:.4f}\".format( endtime*2 ))\n",
    "\n",
    "# 2. Recurrent calculations\n",
    "t = 0\n",
    "dt = 0.0002   #in second, discrete time tick unit, stepping simulator\n",
    "# Vectors for display\n",
    "xx=[]\n",
    "vv=[]\n",
    "aa=[]\n",
    "jj=[]\n",
    "time=[]\n",
    "\n",
    "#Curves\n",
    "while(t < endtime):\n",
    "    t = t + dt\n",
    "    time.append(t)\n",
    "    \n",
    "    if t<t0:\n",
    "        jj.append(jM)\n",
    "        aa.append(jM*t)\n",
    "        vv.append(jM*t*t/2)     \n",
    "        xx.append(jM*t*t*t/6+startX)\n",
    "    elif t<t1:\n",
    "        jj.append(0)\n",
    "        aa.append(aM)\n",
    "        vv.append(aM*(t-t0) + vt0)     \n",
    "        xx.append(aM*(t-t0)*(t-t0)/2 + vt0*(t-t0)+xt0) \n",
    "    elif t<t2:\n",
    "        jj.append(-jM)\n",
    "        aa.append(-jM*(t-t1) + aM)\n",
    "        vv.append(-jM*(t-t1)*(t-t1)/2 + aM*(t-t1) + vt1)     \n",
    "        xx.append(-jM*(t-t1)*(t-t1)*(t-t1)/6 + aM*(t-t1)*(t-t1)/2 + vt1*(t-t1) + xt1)\n",
    "    else:\n",
    "        jj.append(0)\n",
    "        aa.append(0)\n",
    "        vv.append(vM)     \n",
    "        xx.append(vM*(t-t2) + xt2)        \n",
    "        \n",
    "fig = plt.figure(constrained_layout=True)\n",
    "fig.suptitle('Half travel motion S-Curves', fontsize=16, y=1.04)\n",
    "gs = fig.add_gridspec(4, 1)\n",
    "\n",
    "f_ax1 = fig.add_subplot(gs[0, 0])\n",
    "f_ax2 = fig.add_subplot(gs[1, 0])\n",
    "f_ax3 = fig.add_subplot(gs[2, 0])\n",
    "f_ax4 = fig.add_subplot(gs[3, 0])\n",
    "\n",
    "f_ax1.grid(True)\n",
    "f_ax1.plot(time, xx)\n",
    "f_ax1.set_title('x(t)')\n",
    "\n",
    "f_ax2.grid(True)\n",
    "f_ax2.plot(time, vv)\n",
    "f_ax2.set_title('v(t)')\n",
    "\n",
    "f_ax3.grid(True)\n",
    "f_ax3.plot(time, aa)\n",
    "f_ax3.set_title('a(t)') \n",
    "\n",
    "f_ax4.grid(True)\n",
    "f_ax4.plot(time, jj)\n",
    "f_ax4.set_title('j(t)') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}