{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Utilities\n", "\n", "This notebook contains utility code that is useful to multiple notebooks." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Derivative Filter\n", "\n", "In many control schemes that follow trajectories, higher order derivatives are often needed. This section discusses an implementation of a digital derivative based in Laplace system theory.\n", "\n", "Recall from system theory that the Laplace transform is a useful way of analyzing linear time-invariant (LTI) systems. It takes time-domain systems and expresses them in the Laplace domain (a superset of the Fourier frequency domain). The unilateral Laplace transform is defined by\n", "$$\n", "\\begin{equation}\n", "F(s) = \\int_0^\\infty f(t) e^{-st} dt,\n", "\\end{equation}\n", "$$\n", "where $t$ is a time variable and $s$ is the Laplace variable. This transform can be thought as projecting a time-domain function onto exponentials of varying phase.\n", "\n", "This transform is particularly useful when analyzing LTI systems of differential equations. In particular, the Laplace transform of the time-derivative of a continuous system $f(t)$ can be found to be\n", "$$\n", "\\begin{equation}\n", "\\frac{d}{dt} f(t) \\stackrel{\\mathcal{L}}{\\rightleftharpoons} sF(s).\n", "\\end{equation}\n", "$$\n", "\n", "Therefore, differentiation in the time-domain corresponds to multiplication by $s$ in the Laplace domain. Although a pure derivative in the Laplace domain is causal, it is not realizable (i.e., can you build an RLC circuit that realizes $F(s)=s$?). Therefore, a band-limited derivative (i.e., a *dirty derivative*) is used:\n", "$$\n", "\\begin{equation}\\label{eq:dirty-derivative}\n", "G(s) = \\frac{s}{\\tau s + 1}.\n", "\\end{equation}\n", "$$\n", "\n", "From analog filter design, remember that a first-order low-pass filter with cutoff frequency $\\omega_c$ is defined in the Laplace domain as\n", "$$\n", "\\begin{equation}\\label{eq:lpf}\n", "H_{LPF}(s) = \\frac{\\omega_c}{s + \\omega_c} = \\frac{1}{s/\\omega_c + 1},\n", "\\end{equation}\n", "$$\n", "and a first-order high-pass filter is\n", "$$\n", "\\begin{equation}\\label{eq:hpf}\n", "H_{HPF}(s) = \\frac{s}{s + \\omega_c} = \\frac{s/\\omega_c}{s/\\omega_c + 1}.\n", "\\end{equation}\n", "$$\n", "\n", "Note that the cutoff frequency is related to the filter time-constant (and RC circuits) by\n", "$$\n", "\\begin{equation}\n", "\\tau = RC = \\frac{1}{2\\pi f_c} = \\frac{1}{\\omega_c}.\n", "\\end{equation}\n", "$$\n", "\n", "From this discussion, we can see that the dirty derivative \\eqref{eq:dirty-derivative} can be thought of as a filtered version of a pure derivative with bandwidth $1/\\tau$. This is useful as the LPF will prevent the derivative operator from picking up high-frequency components of the input signal and amplifying noise.\n", "\n", "A common value of $\\tau$ is $0.5$ ($20$ Hz). A higher $\\tau$ corresponds to less bandwidth and more rejection of high-frequency components, which results in a smoother output.\n", "\n", "### Digital Implementation\n", "\n", "In order to implement a derivative filter in a computer, it of course must be discretized. The technique used here is to map $G(s)$ to the $z$-domain via the bilinear transform (AKA, Tustin approximation)\n", "$$\n", "\\begin{equation}\n", "s \\mapsto \\frac{2}{T}\\frac{1-z^{-1}}{1+z^{-1}} ,\n", "\\end{equation}\n", "$$\n", "where $T$ is the sample period. Once we have an expression in the $z$-domain, the inverse $z$-transform can be used to give a discrete-time implementation of the dirty derivative.\n", "\n", "#### Resources\n", "\n", "- [SE.DSP: First-Derivative Analog Filter](https://dsp.stackexchange.com/questions/41109/first-derivative-analog-filter)\n", "- [Blog: Causal, but not Realizable](http://blog.jafma.net/2015/10/04/differentiation-derivative-is-causal-but-not-exactly-realizable/)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "class DirtyDerivative:\n", " \"\"\"Dirty Derivative\n", " \n", " Provides a first-order derivative of a signal.\n", " \n", " This class creates a filtered derivative based on a\n", " band-limited low-pass filter with transfer function:\n", " \n", " G(s) = s/(tau*s + 1)\n", " \n", " This is done because a pure differentiator (D(s) = s)\n", " is not realizable. \n", " \"\"\"\n", " def __init__(self, order=1, tau=0.05):\n", " # time constant of dirty-derivative filter.\n", " # Higher leads to increased smoothing.\n", " self.tau = tau\n", " \n", " # Although this class only provides a first-order\n", " # derivative, we use this parameter to know how\n", " # many measurements to ignore so that the incoming\n", " # data is smooth and stable. Otherwise, the filter\n", " # would be hit with a step function, causing\n", " # downstream dirty derivatives to be hit with very\n", " # large step functions.\n", " self.order = order\n", " \n", " # internal memory for lagged signal value\n", " self.x_d1 = None\n", " \n", " # Current value of derivative\n", " self.dxdt = None\n", " \n", " def update(self, x, Ts):\n", " # Make sure to store the first `order` measurements,\n", " # but don't use them until we have seen enough\n", " # measurements to produce a stable output\n", " if self.order > 0:\n", " self.order -= 1\n", " self.x_d1 = x\n", " return np.zeros(x.shape) \n", " \n", " # Calculate digital derivative constants\n", " a1 = (2*self.tau - Ts)/(2*self.tau + Ts)\n", " a2 = 2/(2*self.tau + Ts)\n", " \n", " if self.dxdt is None:\n", " self.dxdt = np.zeros(x.shape)\n", " \n", " # calculate derivative\n", " self.dxdt = a1*self.dxdt + a2*(x - self.x_d1)\n", " \n", " # store value for next time\n", " self.x_d1 = x\n", " \n", " return self.dxdt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Rotation Matrices\n", "\n", "Throughout this notebook, we use *right-handed*, *passive* rotations to express the same geometrical vector in various coordinate frames. When Euler angles are used, the *intrinsic* 3-2-1 (Z-Y-X) sequence is used. This means that a yaw ($\\psi$) rotation is first performed in $\\mathcal{F}_A$ around the $\\mathbf{k}^A$ axis to get to $\\mathcal{F}_B$. Then, a pitch ($\\theta$) rotation is performed about the $\\mathbf{j}^B$ axis of frame B to get to $\\mathcal{F}_C$. Lastly, a roll ($\\phi$) rotation is performed about the $\\mathbf{i}^C$ axis to get to $\\mathcal{F}_D$. The order of rotation is important, and if data from $\\mathcal{F}_A$ is meant to be expressed in $\\mathcal{F}_D$, the rotations are composed as\n", "$$\n", "\\begin{equation}\n", "x^D = R_C^D(\\phi) R_B^C(\\theta) R_A^B(\\psi)x^A = R_A^D(\\phi,\\theta,\\psi)x^A = R_A^D x^A.\n", "\\end{equation}\n", "$$\n", "\n", "Python implementations of these passive rotations are defined below." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# input angle in radians\n", "def rotx(ph): return np.array([[1,0,0],[0,np.cos(ph),np.sin(ph)],[0,-np.sin(ph),np.cos(ph)]])\n", "def roty(th): return np.array([[np.cos(th),0,-np.sin(th)],[0,1,0],[np.sin(th),0,np.cos(th)]])\n", "def rotz(ps): return np.array([[np.cos(ps),np.sin(ps),0],[-np.sin(ps),np.cos(ps),0],[0,0,1]])\n", "def rot3(ph,th,ps): return rotx(ph).dot(roty(th).dot(rotz(ps)))\n", "# input angle in degrees\n", "def rotxd(ph): return rotx(np.radians(ph))\n", "def rotyd(th): return roty(np.radians(th))\n", "def rotzd(ps): return rotz(np.radians(ps))\n", "def rot3d(ph,th,ps): return rot3(np.radians(ph),np.radians(th),np.radians(ps))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## PID Controller\n", "\n", "One of the simplest forms of feedback control is the proportional-integral-derivative (PID) controller. It is perhaps the most widely used form of control because it is an intuitive technique for minimizing error. However, stability can only be guaranteed for second-order systems." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Digital Implementation\n", "\n", "In order to implement a PID controller on a computer, we need to discretize our continuous expression for the integral and derivative terms of the PID controller." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "class SimplePID:\n", " def __init__(self, kp, ki, kd, min, max, tau=0.05):\n", " self.kp = kp\n", " self.ki = ki\n", " self.kd = kd\n", " self.min = min\n", " self.max = max\n", " self.tau = tau\n", "\n", " self.derivative = 0.0\n", " self.integral = 0.0\n", "\n", " self.last_error = 0.0\n", " \n", " @staticmethod\n", " def _clamp(v, limit):\n", " return v if np.abs(v) < limit else limit*np.sign(v)\n", "\n", " def run(self, error, dt, derivative=None, pclamp=None):\n", "\n", " # P term\n", " if self.kp:\n", " # Proportional error clamp, it specified\n", " e = error if pclamp is None else self._clamp(error, pclamp)\n", " p_term = self.kp * e\n", " else:\n", " p_term = 0.0\n", "\n", " # D term\n", " if self.kd:\n", " if derivative:\n", " self.derivative = derivative\n", " elif dt > 0.0001:\n", " self.derivative = (2.0*self.tau - dt)/(2.0*self.tau + dt)*self.derivative + 2.0/(2.0*self.tau + dt)*(error - self.last_error)\n", " else:\n", " self.derivative = 0.0\n", " d_term = self.kd * self.derivative\n", " else:\n", " d_term = 0.0\n", "\n", " # I term\n", " if self.ki:\n", " self.integral += (dt/2.0) * (error + self.last_error)\n", " i_term = self.ki * self.integral\n", " else:\n", " i_term = 0.0\n", "\n", " # combine\n", " u = p_term + d_term + i_term\n", "\n", " # saturate\n", " if u < self.min:\n", " u_sat = self.min\n", " elif u > self.max:\n", " u_sat = self.max\n", " else:\n", " u_sat = u\n", "\n", " # integrator anti-windup\n", " if self.ki:\n", " if abs(p_term + d_term) > abs(u_sat):\n", " # PD is already saturating, so set integrator to 0 but don't let it run backwards\n", " self.integral = 0\n", " else:\n", " # otherwise only let integral term at most take us just up to saturation\n", " self.integral = (u_sat - p_term - d_term) / self.ki\n", "\n", " # bookkeeping\n", " self.last_error = error\n", "\n", " return u_sat" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Numerical Integration\n", "\n", "Before our discussion of the dynamical equations of motion that describe the time evolution of the quadrotor in $\\text{SE}(3)$, we should discuss how to properly integrate these differential equations numerically for use in a simulation. The two most common types of numerical integrators are Euler's method and the fourth order Runge-Kutta method, or RK4.\n", "\n", "### The Euler Method\n", "\n", "The Euler method is an explicit algorithm that uses the limit definition of the derivative:\n", "$$\n", "\\frac{df}{dt} = \\lim_{h\\to0} \\frac{f(t + h) - f(t)}{h}.\n", "$$\n", "Let $g(t) \\triangleq \\frac{df}{dt}(t)$. Disregarding the limit and rearranging the difference quotient yields\n", "$$\n", "f(t+h) = f(t) + g(t)h.\n", "$$\n", "Using the step size $h$ as the sample period and using discrete notation, we have\n", "$$\n", "f[k+1] = f[k] + g[k] T_s.\n", "$$\n", "\n", "### RK4\n", "\n", "Of course, the Euler Method is a first-order, rough approximation of differential equations. Given a differential equation $\\dot{y} = f(t,y)$, the RK4 method is given by\n", "$$\n", "\\begin{align}\n", "y_{n+1}\t&=\ty_{n}+\\frac{h}{6}\\left(k_{1}+2k_{2}+2k_{3}+k4\\right) \\\\\n", "t_{n+1}\t&=\tt_{n}+h,\n", " \\end{align}\n", "$$\n", "\n", "where\n", "$$\n", "\\begin{align}\n", "k_{1} &= f\\left(t_{n},y_{n}\\right) \\\\\n", "k_{2} &= f\\left(t_{n}+\\frac{h}{2},y_{n}+\\frac{h}{2}k_{1}\\right) \\\\\n", "k_{3} &= f\\left(t_{n}+\\frac{h}{2},y_{n}+\\frac{h}{2}k_{2}\\right) \\\\\n", "k_{4} &= f\\left(t_{n}+h,y_{n}+hk_{3}\\right).\n", "\\end{align}\n", "$$\n", "\n", "For an expository derivation of RK4, see [2]." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def rk4(f, y, dt):\n", " \"\"\"Runge-Kutta 4th Order\n", " \n", " Solves an autonomous (time-invariant) differential equation of the form dy/dt = f(y).\n", " \"\"\"\n", " k1 = f(y)\n", " k2 = f(y + dt/2*k1)\n", " k3 = f(y + dt/2*k2)\n", " k4 = f(y + dt *k3)\n", " return y + (dt/6)*(k1 + 2*k2 + 2*k3 + k4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### RK4 Example\n", "\n", "Consider the autonomous evolution of position $x \\in \\mathbb{R}$ as\n", "$$\n", "\\dot{x} = \\cos(\\omega(t))\n", "$$\n", "\n", "where\n", "$$\n", "\\omega(t) =\n", "\\begin{cases}\n", " 2 \\pi t, & 0 < t < 2.5 \\\\\n", " 4 \\pi t, & 2.5 \\leq t < 5 \\\\\n", "\\end{cases}.\n", "$$\n", "\n", "Note that this equation is not smooth (there is a jump at $t = 2.5$), but we are still able to numerically integrate it." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsvXl0W/d17/vdAAiAAAFwBkeJkiXLlhRHjumhqePQGWqnyYt9V4amt7cv6Uuf27eavrRZbZrcNmmbtKvp7b1N3n23XY1v2zQd7kvSDLXbOnEcW/SQeJJkS5ZkTZYoiRPACRNBgBj2++OcA4IUSJwZONTvsxYWgYNz8Pt9ccDf/g177x8xMwQCgUAgMIqr0RUQCAQCwfZAGBSBQCAQmIIwKAKBQCAwBWFQBAKBQGAKwqAIBAKBwBSEQREIBAKBKQiDIhAIBAJTEAZFIBAIBKYgDIpAIBAITMHT6ArYSXd3N4+MjOi6dnl5GcFg0NwKNTlC8/WB0Hx9YETz0aNH55m5p95515VBGRkZwZEjR3RdOz4+jrGxMXMr1OQIzdcHQvP1gRHNRHRZzXliyksgEAgEpiAMikAgEAhMQRgUgUAgEJiCMCgCgUAgMAVhUAQCgUBgCsKgCAQCgcAUhEERCAQCgSkIg2Ix52NpfOVH5/D8GwuNroptnLuONb9w8frRLBBs5LoKbLSbE5MJfOivnke+WAZwHl/+uTfjP9w61OhqWcpGzV/5uUN48NbBRlfLUo5fTeBDX30eq7Lm/+cjh/DAoe2tWSCohRihWES5zPj0t0+gK+jFs5++F3fs6sTn/+UUlpZXG101y7hG80gnPvcvJ5HIbl/NJVlzT5sPz376Xtw+0oHf2+aaBYLNEAbFIp54PYYzs2n8zntuwnBnAH/04EFkVov42k8mGl01y/jh6VmcmU3jMz97M4Y7A/iiovnHE42ummX88NQszsbS+Ix8n7/44EFk8kV8/SeqMlUIBNsKYVAs4hsvXUFf2I/3vqkfAHBjNIR79vbgWy9fRanMDa6dNfx/L11Ff2RN876+EN62twffOrKNNb98FYPtrfhZWfNNfWHcvacb33z5yrbVLBBshjAoFpDKFfDM+Xk8eOsgPO61r/hDo0OYTeXw8sRiA2tnDclsAc9dkDS7XVQ5/qHbhjCTzOHIdtV8fg4P3jqwXvPoMKaTORy9vNTA2gkE9tNQg0JE9xPRWSK6QESfqfH+PUR0jIiKRPTBDe+ViOhV+fGofbWuz3Pn51EqM951c++64/fu64XX48Ljp2YbVDPrePbCXG3NNymaYw2qmXU8c34OZQbeeXN03fF33NQLr3t73meBYCsa5uVFRG4AfwHg3QAmAbxMRI8y8+mq064A+BiA36rxESvMfMjyiurgqTNxRFpbcGi4fd3xoM+Dt97QhWfPzzeoZtbx1Jk42gMtODTcse54m8+Du3Z34dnzcw2qmXUcPhNHR6AFbx5af5/bfB7cubtzW2oWCLaikSOUOwBcYOaLzLwK4BsAHqg+gZknmPkEgHIjKqiHcpkxfnYO99zYs266S+HOXV24EM9gPpNvQO2soVxmPH12Dvfs7Vk39aNw1+5OnI9nsLDNNI+fm8Pbb9xMcxfOxTJY3MZefQLBRhppUAYBXK16PSkfU4ufiI4Q0QtE9KC5VdPPedlYvP3G2pub3bW7EwDw0qXts6ZwNpbGwvLqpprv3NUFYHtpPjObxuLyKt6+r959FoGOgusHJwc27mTmKSLaDeApInqNmd/YeBIRPQTgIQCIRqMYHx/XVVgmk1F17dNXCwCA4uw5jI9fuOb9YpnhcwPfefYEAgtnddXFLtRqHleh2esGvv3sCbRuE82Hr0iaCzPnMJ7cQvMzJ+Cf3x6atxNCszU00qBMARiuej0kH1MFM0/Jfy8S0TiAWwFcY1CY+WEADwPA6Ogo690CU+32mT/4zgm0B2bxcz97L4iunQoBgDsuvYjJVB5jY/foqotdqNX82LePoyMQw4e30HznpRcxld4+mv/9n4+jMxjHh96zxX2++CImM9tH83ZCaLaGRk55vQxgLxHtIiIvgI8AUOWtRUQdROSTn3cD+GkAp7e+yh5euZLAoeH2TRsZALhtZwfOxdNYzhdtrJl1vHq1vua37OjA2Vga2dXtofmVqwncWk/zzg6ci6WxslqysWYCQeNomEFh5iKATwB4HMDrAL7FzKeI6AtE9H4AIKLbiWgSwIcAfJWITsmX3wzgCBEdB3AYwJc2eIc1hEy+iHPx9DXeXRs5OBABM3BmNmVTzawjnSvgfDxzjXfXRg4OSppfn0nbVDPrSK4UcCGeqXufDwyEUWbg9W1wnwUCNTR0DYWZHwPw2IZjn696/jKkqbCN1/0EwJssr6BGTlxNgBn1G5rBMADg5FQKt+3stKNqlnFiMilp3lG/cQWAU9NJ3LZza+PT7JyYTACor/ngYAQAcGoqibfscLZmgUANIlLeRE7PSD3RN8kNyWb0hf3oCnpxajppR7Us5fS0Os39ET86g16cmnJ+b12t5oGIH+2BFpyadr5mgUANwqCYyJnZNLrbfOhq8215HhFh/0AYJ7dB43pmNo2ekA+dQe+W5xERDgyEcXIbGNGzs2lEwz60B64fzQKBGoRBMZGzs2nc1BdSde6BgQjOx9PyHhrO5WwspVrz/oEwzsWcr/nMbBr7+sKqzj04EMG52QwKJWdrFgjUIAyKSZTKjPPxNPZpaFwLJcaFeMbimllHqcw4H8tgX1Sl5n5J8xtzztVcLJVxYS6jyYiulsqO1iwQqEUYFJO4sphFrlBW3bju7W0DAFxwcEMzsbCMfLGs2oju7ZXOc7IRnVjIYrVYxo0q7/MNPfJ9drBmgUAtwqCYxFnZNVRt47qrOwgXObuhOTsruQDfpHL6Z3dPELRtNKs3KE7XLBCoRRgUkzg7mwERVPdc/S1uDHcG8IaDG5qzs2m4CNgbbVN1vr/FjeGOgKNHZWdjkuY9veo0t3rdGGxvFQZFcF0gDIpJvDGXwWB7K1q9btXX7O1tc3RD88ZcBkMdAfhb1Gve09vmaCP6xlwGw53aNDv9PgsEahEGxSQuzS9jd4+6XqvCDb1tuDS/jKJDPYAkzUFN1+zpbcPF+WXHbo97aW4Zu7uvL80CgVqEQTEBZpYaV60NTU8bVktlXF1asahm1qFo3qVHc7GMq4tZi2pmHWuatXUc9vRKmieXnKdZINCCMCgmMJfJI5Mvam9ce53rARRP55FdLWk2ojc4WHMslcdKoYRdOkZlgDM1CwRaEAbFBC7NLQOAZoOyW+7pTswvm14nq7lY0axxmk9ujCcWHKh5XjIIWo2ocp8vOfA+CwRaEAbFBJSGQqtBiQRa0B5oweVF5zU0Fc0ae+vtAS/Cfg8uLzhv+kfvfW4PtCDk9+CKA6f5BAItCINiApfml+H1uDDQ3qr52p2dAYc2rhn4PC70h/2arx3pDjpyhHJpbhn+Fhf6NGomIox0BTHhwPssEGhBGBQTuDi/jJGuANyuzTdb2oydXQ5tXOUFeZcOzTsca0SXMdKlU3NXAFcceJ8FAi0Ig2IClxekhkYPO7sCmFpacVzCxImFLHZ2BXRdO9IVxFRixXEJEycWtHu1KYx0BTC5tOJYF3GBQA3CoBiEmXF1cQXDnfoa151dQZQZmEo4x3WYmTG5lMVwh17NAZTKjCkHuUtLmo3d52KZMZ3ImVwzgaB5EAbFIAvLq1gplDDUoX39BECll3/ZQdMhc5k8coWyAc1SL/+ygxap59J55IsGNMuGyInTmwKBWoRBMcik3Ms20lsH4Kg1hYpmnb31EQca0asG7/NIt/OMqECgFWFQDKJEfOttXHvafAh43Y7quRrWHPKhtcWNiXnnNK5KlLveEUpvyAd/iwuXRSyKYBsjDIpBrhpsaIgIg+2tjlpPUEYogzrcpAFZc0crphLOMSiKER3SOUKp3GcHrZUJBFoRBsUgk0sr6Ax6EfR5dH/GUIezGprJpSy6rjvNK+hu82nKJr2RoY6AozQLBFppqEEhovuJ6CwRXSCiz9R4/x4iOkZERSL64Ib3PkpE5+XHR+2r9XquLmYxrHN0ojDY0Vrp9TuBq4srGNI53aUw2O4wzUtZ3aNQBafdZ4FAKw0zKETkBvAXAN4DYD+Anyei/RtOuwLgYwD+14ZrOwH8PoA7AdwB4PeJqMPqOtdicmlF9zSIwmB7AMmVAjL5okm1shazGtdEtoBlp2g24BquMNjeisXlVWRXnaFZINBKI0codwC4wMwXmXkVwDcAPFB9AjNPMPMJABujwe4D8AQzLzLzEoAnANxvR6WrKcuxFEOdxhtXAI5YRymVGdOJFd3eTgrK+osTpoAUzUaNqHL9tAM0CwR6aKRBGQRwter1pHzM6mtNI57OY7VUNrFxbf5F6lgqh0KJMWzQiA45yIjOpnIoltm0+yymvQTbFf2rqg6BiB4C8BAARKNRjI+P6/qcTCZzzbXnlkoAgMWr5zGeu6S7jks5aQB2+KUTcM226P4cs6ml+eyirPnKeYyvmKH5OKjJNZ9ZrLrPKxd1f/aiovnF48BMc2ve7gjN1tBIgzIFYLjq9ZB8TO21YxuuHa91IjM/DOBhABgdHeWxsbFap9VlfHwcG69demUSePE43jt2F27QuP1vNeUy43ee/QECPUMYG7tZ9+eYTS3NC0cngZckzVq3PK6mXGZ8+tnvI9AzjLGxmwzW1DxqaZ6v0qw3lxcgTZ19+pnvI9jb/Jq3O0KzNTRyyutlAHuJaBcReQF8BMCjKq99HMDPEFGHvBj/M/IxW7m6aCweQ8HlIvS3+x0x/aPE3QwaXE9wuQj9EWe4Dl9dzIIIGGjXnqq/GrdynzVqzhVKeObcXEOSaZ6ZTTVkYzBFcyOSaV6Pms2iYSMUZi4S0ScgGQI3gL9l5lNE9AUAR5j5USK6HcD3AHQA+N+I6A+Z+QAzLxLRFyEZJQD4AjMv2q1hcimLnpAP/hb9sQkKTgl6m1paQW/IB5/HJM0O2Gd9KmG2ZvX3uVRm/NxXn8fxySTuubEHX/+l20GkPX2+Hh4/NYtf/cej8LgI//TLd+GOXZ22lFsslfHhrz6PE5NJjO3rwdc+Zp/mH5ycwf/1T8caovlDf/U8XptK4t59PfhbGzWbSUPjUJj5MWa+kZlvYOY/lo99npkflZ+/zMxDzBxk5i5mPlB17d8y8x758bVG1H8mmdO1qVYtnBKXMZvKoT9irKeu4JS4jFgqh76IWfc5oEnzo8encHwyiX3REJ45N4dnzs+bUo96lMuMP3nsdTk1kAdf+v7rtpQLAI8en8aJySRujLZh/OwcnrVJc6nM+JPvn0FvSNL8pz84Y0u5APDIq9N4bUrSfPjsHH58YcG2ss1ERMobYDqxggGTGtehjgDm0nnkCiVTPs8q4qk8ojp2aazFUEcr4uk88sXm1hxL5dAX9pnyWUMdrYilc6r3v/nusSmMdAXwr79+NzoCLfjesUlT6lGPV64uYWIhi//8szfjE/fuwbErCduSeX732BR2dQfxr79+N9oDLfjeK2qXVo1x7MoSLsuaf+3eG3D08hKu2JS09buvTGK3rDnS2oLvvmLPfTYbYVB0wsyYSebQb1bPVV6TmEk2934Zs6mcaQZFWXuaafI9QmaTJmruaAWz9Jn1SOUKeP6NBdx3sA9ejwvvvDmKJ8/EbVlL+eGpGFrchHfc3Iv7D/YBkKbArCa5UsALFxdw34E++DxuvPOmKJ58PWbLusITp2Pwul14x029eM/BfgDAD0/boDlbwAsXF3H/QUVzL558Pe7ItRRhUHSSyhWRXS0ZXqhVqMSiNPEUUK5QQnKlgKhJvfVKQGcTrx2trJaQyhXNG5UpsSgqYo6OTiyhWGa8/cYeAMDbb+xBOlfE6zMpU+qyFc9fXMCtOzoQ9rdguDOA3T1BvHDR+mXKo5cX12m+58ZupHJFvD6Ttrzs599YwFt2tiOkaO4O4oWL1k89vTyxiFKZcY9yn/f1ILlSwJlZ6zWbjTAoOplJSo1gn2lTXs0f3BhP5QHAxMZVChRsZiMaT0sjCTNHKIA6zUcuL8LtIhwabgcAjI5I2YWOTCyZUpfNyK4WcWo6hdtH1rIZje7swNHLSyiX2dKyj0wswVOl+fYRaVH8yGVrjdlyvojTMymM7lxbhL9N1sxssebLS2hxE9481F4pFwCOTNjuZ2QYYVB0okzTmDXlpTRYzTzlNZsyt3GNRqSRTlNrTiqazRmVKd+dmimvV64ksL8/jIBXcsbsj7RisL0Vx65Ya1Bem0yiVOZKwwYAozs7kVwp4KLF7rSvXElg/0C4ktV5oL0VAxE/jl1JWFrua1M1NI90YClbsNyF+JUrS9g/EKloHmxvRV/Yes1WIAyKTqblEYpZU15ejwvdbV7EUs3buMZMNig+jxtdQW/FUDUjsbS5ozJ/ixudKjQzM07PpHBwMLzu+P6BsOVTXqflzz84EFlXLgBLy1Y0H6gqVyo7Yr3maenzD1R930o9rJxuq9zngbVyiQgHbLjPViAMik5mkzm4SNpx0Sz6Iv6m7q2vGRRzNc8mm3jKS9EcMsegAEBf2F93hDKbyiGRLeDm/vUG5eb+MC7NL1vqDfj6TApdQS96Qmv3eW+0DR4XWdrIzSRzSK4UsL8/tO74/v4QLs5lLNfc3eZFb9V93tPbBrfFmqcSK0jnijXv80WL77MVCIOik+mE5PnjcZv3FappaBpJPJ2Hz+NCpNW8PFR9YT9m5bWZZiSWysHf4kK41bwY4L6Iv+4IRWnENjY0+/tDKDNw1sIF2zOzadzUH1oXWOfzuHFDT5uljetmmm/qD6PMwLmYdZpfn01dU66/xY0beoKWaj4jj35u3mBEb+4Po1RmnI9lLCvbCoRB0clMcsW0AD8FNQ1NI1HcZ82M4I02+QhlVo67MVWzio7DxTlp3n5v7/p8aXt6pYbn4rw1DQ0z4+LcMvb2hq55b0+0zdL1hDXN68tWvgOrylY07+m9Njfd3t6QpetGyn3cs1FztG3d+05BGBSdmBmDotAX9iORLTTtMDeWypk63QUA/WE/lppds4nTXQDQH/FjYXl1y4DOS/PLiLS2oD3gXXd8uLMVLgIuzVvjDTifWUUmX8RI17Wp+nd1BXF1acWyOJhLC8toD7QgElg/Ah7uDIDIOoMyl84ju1rCSNe1iT93dgVwdTFrWUzIpfksOoPea0b9O+TN3CYsus9WIQyKDqSgRitGKJKBatZpr3javCh5haj8HcabdNornspV6mgWfeH6micWljFSI7Oxz+PGQHsrJixqXCfkaPidNcre2RVAqcyWpcuZmF+u2aj7W9wYiFinWTFUtb7vke4gihZr3lnDeEua/ZX74RSEQdGBNIooo9+kPF4KSkPTjNNezGxqxLiCYpRnmnDai5mlzAAhc0dlSuzSVvd5Yj6LXTUaGgDY1R20rKFRGtddNRp2JXW/ZcZsfnnT7QFGugO4ZFEaFOW73ErzJYu+74mF5ZrlApIxa0TWYyMIg6KDisuwBWsoQHOOUNL5IlYKJdOnvJrZiKZyReQKZdONaF9k65ijXKGE6eRKzR4zAIx0SQ2NFQF3E/PL8Lio5nbHSn2sMGaS5lzNEQogabZuVJZFi5tqhgAo9blsQdkrqyXMJHOb3+fuoG3508xCGBQdKEGNZkXJK6jpuTaKuMkxKArNbEQVzb1mG1FZc2wTzVcXs2DGpo3rzq4A0rkilrIFU+sFSMZiuDNQ03uxK+hFyOexpGG/siiNPka6Nx+VJVcKWFpeNb3sifllDHfU1tzd5kXQ68aEBaOjy4ubT7UBwEhXAEvZApIW3GerEAZFBzNyQ2NW6nqFNp8HIZ+nKRvXmMlpVxRC/hYEve6mNKKK5j6zNfs8CHjdm45QtprTB6qmYSxo2C/NZ2suyANSwN1Id9CSqaeK5i1GKIA1U0+X5muvVwFVmi34rie2mF4ErNVsFcKg6CCekoIau00MalSQ3Gibr3FdS0FibuMKNLFmi0ZlRIS+sH/TrAhbzekDVVNPJjdyzIzLmzgDVJdtxQhloo4RtVZzdlNDppRtxTSf4qm31agMsG7NygqEQdFBLJVDT8gHt8v8HdWkQL/ma1xjafOj5BX6mzT+xuxUM9VsFXN0ZTGLSOu17rMKwx2SG+3lRXNHCvOZVWRXS9jZWbuBA6RpmKnEiulutFcWs5LL8CZBs8OdrZJmk0dHc5k8Vgqlmp5WCiNd0qZoJZMTY15ZlFyGQ/7NNEt1MluzlQiDooOYiZtMbaSvSXvr8VQeIb+nkqjQTKJh/6brCY0knsoh7PdUkvaZyVZZEaYTucp2BrXwelzoDfkwY3La/2n58wY7Nm9cB9pbUSoz4mlz3bynEytbavZ53Ohp85nuDTgtr4duVfaaZnN/o/U0+1vc6LZAs5UIg6KDWCq3LuePmfSF/ZjL5Jtucx0rXIYV+iN+xNJ503uARjFzM7GN9EWkKa9a6eCnEyt11+cG2lsr3oZmoRiUrRKeKvWaNt2Y1d9Oe6C9tWIAzCtX0by1Qak+18yy6yWXHWz3N/V+QRsRBkUH8XTedM8fhb6IH6UyYz5jvjeLEWJp86PkFfrCkuaFTHMFN1o9Ei2WGfPL12qeSqxgsE5DY0XjqjRcW/WalXqZ3cjV660r9bKiUVc+e6tyAWDKxO+bmdV3HIRB2b6sFstYXF41PR2HQrPGZcRTees0R5pz++N4KmddxyGsuA6vNyipXAHpXLFuQzPY3oqpxIqpsSgzyRwCXveWyT+VdENmGrNUroB0vli3tz7Q7sd00lzNU4kVBL3uLZN/KsG3ZjbsqVwRy6ulukZU6ThYvcmXWQiDopG5jOI+a90IBUBTJUwsy/PHZqcgUWhGI1qW1wnMdhlW6NskQ4AS41S35xrxY7VYxoKJcRlKj3mrRJhBnweR1hZTG1c1007K+7lC2dT4GzWaQ/4WhP2ehmleKZSQcEgsSkMNChHdT0RniegCEX2mxvs+Ivqm/P6LRDQiHx8hohUielV+/JVddbbS8wdozkC/xewqCiU2PQWJQiXQr4kMysLyKopltvw+b9SsduO2fgvm9acT6vLTmT0No3b307XRkYllJ3OqUiiZrVn5rHrft5KNwynrKA0zKETkBvAXAN4DYD+Anyei/RtO+ziAJWbeA+DLAP606r03mPmQ/PhVWyoN66KnFToDXnjdrkrwZDNgtRHtCnrR4qammvKyYjOxarqDPnhc12pW23MdbDd/6mmqjnfZWtl+TJt4r9Ss3VS/b2bjOq1ivQowf81KzdqNUi7QfNPBm9HIEcodAC4w80VmXgXwDQAPbDjnAQBfl59/G8A7ycyNKXSgRE9b5eXlchF6w76mcqNVsuJaNeXlchF6Q83lOhxPW2tEJc2+a6b5phMrcMvfx1aY7XmUL5Ywn8mryv7QHzG/t+5x0bodImuhjNrMcpfOFUqYz6xiQMU2FMr6jVlMJ3NocVPd4GirPMysopEGZRDA1arXk/KxmucwcxFAEkCX/N4uInqFiJ4mordZXVmFeDoHt4vQFfTWP1kn/U22FbBVEePVNNv2x7NJa1LNVFMr5mg6kUNf2F83aLYj0AJ/i8u0hkaphxqDMtDeiuRKAZl80ZSypxMr6IvU19wZ9MLncZk2OprRqDmRLWDZRM39kVa46mjuCnrh9Zh3n62mbpQaEf0UgP8E4G0A+gGsADgJ4N8B/CMzJy2tYW1mAOxg5gUiug3AvxDRAWa+Zq9OInoIwEMAEI1GMT4+rqvATCaD8fFxHD+XR8QLPPPM0waqvzWUy2EiVdZdV7NQNL94QVr4ff3YCzhvQXYAAHDlc7jUTJrPS5pPH3sB5yzUfHFhvebTEysIElR9D+1exqvnr2B8PG6oHplMBo+NvwAAiE+cxXj6wpbnJ6elRvWRHz6DwZDxPunpyysIQJ3mDi/jlXOXMR6IGSozk8ng+09r1/wvTzyDwTbjml+/vIJW2K/Z6v+vLQ0KEX0fwDSARwD8MYA4AD+AGwHcC+ARIvpzZn5UR9lTAIarXg/Jx2qdM0lEHgARAAss+dDlAYCZjxLRG3KdjmwshJkfBvAwAIyOjvLY2JiOqko3fmxsDH/zxosYdhUxNvbTuj5HDc9lTuO1F6/g7W9/u6lbz2pF0fz44mvoCs7iXe+417KynkmfxumXr0Dv/TGLNc0n0B2LWap5PHUKZ49OrtP8+y8fxi1D7Rgbu7Xu9bvPv4BcoWT4tzg+Po6+wT3Ay8dx/z13YnfPtVvhVtN6cQF/deIF7LjpTXjb3h5DZQPA5186jEPD6jTvOvcCCqUyxsbeaqjM8fFxRPtlzW+/a9N9WBT8Fxfw1RMvYOe+W3D33m5DZQPA5156Crft6FCleeTs8ygzm6LZ6v+veqb2F5n548z8KDNPM3ORmTPMfIyZ/xszjwH4ic6yXwawl4h2EZEXwEcAbDRMjwL4qPz8gwCeYmYmoh55UR9EtBvAXgAXddZDE/FUHr0WeTspRMN+rBRKSJs0vDZK3MKIcYVo2Ifl1ZJp0yhGsTKoUSEa9iOdLyK7KmlmZjneR93vS0owaU4wqLJO1qtC85qHmvGymWWXdJXOD30RfyWvnFGU9DFq/p8rcUMmOMtU7rPK35eUVaG5gn43Y0uDwszz1a+JKExEncqj1jlqkddEPgHgcQCvA/gWM58ioi8Q0fvl0/4GQBcRXQDwKQCKa/E9AE4Q0auQFut/lZkX9dRDK1p+/HpRPMiaZZF61oK95Dei/HM1i7u0lalmFJSGTGnM1zYxU1dub9iPeNqcoLdYKo82nwdtvvq52hSHATMaV62bmPWGfYil8iZpziHk8yCoRrPyP2mCMUutFJEvllUZb0DOdZdyRnCjqkx/RPQrAP4QQA6AoooB7DZSODM/BuCxDcc+X/U8B+BDNa77DoDvGClbD/liCUvZgmUR4wprvaE89kZDlpalhlgqjzcNRiwtQ2lQ4qkc9vRuPeViB/F0Dm8etkdzLCXt2qfVJT0a9qFQYixlC+g06CQSS6vPCtDqdSPs91Tqa4Q1zSob15AU0JlcKaA9YExzXIPmgNeDkN9TMf5G0Jq5uzfkQ756XVxKAAAgAElEQVRYRmqluGkG6mZBberY3wJwUO9oZLsQt2iTqY1ETRxeG6VQKmNh2Y7pH/N6gEYplMqYz6xaP0KpaJZ+V1o3Mav+nRg1KPFUTlNHKWrSdFtFs8ppvmhVZ8uoQdE6ramMFIyiNa6rojmda3qDotZd4Q0AzknKbxFKbIJVQY0Kyuc3QyqS+UwezPYZ0dlk4+eK59I2dRxCa6MyYK2hUbtGVzHCpjRy2hKeRsPmrGVob1zN1JzTtB4aDftM+64BLfe5eTqY9VA7QvksgJ8Q0YuQvasAgJn/b0tq1aRYHdSosDa8bvwPaG2nRmuNaFDe/rgZ/mlmLY6SVwi3euD1uCqLw5VFYrXrCRWDZMwIMzNiGh0vesM+vHjR+E6CMY2dNLMaV8kZQOMIJeTHi5eML9WudRy0GtHGd7bqodagfBXAUwBeA9BcG3XYSNymhgYw14PHCFbtJV+L3rDP9E2M9BDX+A+vFyJCNOxbN0JRuzAOVE2ZGWxcs0VIi8Saeutr+7nUC87bingqj5BP/cZtSjS9Uc3LBSlzuFrjDax3gjDizq914zYznSCsRq1BaWHmT1laEwcQS+fR4iZ0GJy7VYNZUwpGsToFSTVmzcsbRalDn0WpZqrpDa1pjmucdvJ53OgItBj+nSTykp+Ntt66D8UyYzG7Wjd9yFZoWRgHpF0M2wMthn8na5q1TXmZ4QShdWSkOEE4waCoXUP5PhE9RET9G92GryeUnRqN9MjU0iz5vGaT1qeaUYhusS2uncymcvC4CJ22dBzWRmUxjQvj0vXGjXAip92gmJUhWk+8TzRkfHE8kZcmWrQuygNmaNbukq7s8NnsqDUoPw95HQXAUflxTVT6dkdrD9IIfWE/4ul8zS1i7SQmB3LaYUSjJsZVGEFZrLWl4xDyV9ZAtLjuVq4P+w2vtSmNq5YpL2WqyOj6jZ7GtTfsq3jG6WVJNqJaF+UBc4yo1uDoZhm910OVQWHmXTUehmJQnIhWrxAjRMP+ypRCI5ECOa2f+gHWTyk0kngqb1lm5Y30hn2VaHl9vXWf4YZmSZ7+0erlBRhrXJWIca1GNGqKEVUMioY1FBOcIJTMAFrWbpSym8FJpx5bGhQiurvO+2EiOmhulZoXrXOfRlB6Q42eAorZECWv0CzR8nqmnvSiNFLnYxlpkVhHz3Uuk0fJwEg2kWOE/OoXxgGgp82451EiW8Bqqaxjms9nePSeyLOmhXHAHCeIpWxB2qxOsxE1rtkO6o1QPkBEPyGizxPRe4noDiK6h4j+DyL6BwD/BqB+7udtwGqJkVwp2GZQKlMKDV6YtyMFiUJ1AFcjsSPVjIJSzmtTSfm19sa1VGYsLBto2PPad6b0elzoCnoN3auYToePaNgva9Y/etej2QwnCL2b1TXLjEU9tuySMPNvyovvH4CUAkVJX/86gK8y83PWV7E5WBsi27eGAjTW9zxfYqRyRdtHZY0c2ueLjHSuqHlKQi/KCOW1SX0GpXotQ6+bcyLPiHZr/10bXb+pxHVpXTeqcqOttynXZiRyjP4e7d+X0bUMvTuBVq/fGPGqs5q6Y1w56eL/lB/XLRWDYlNDY5a/vRGSOtxJjbCmuXFGVLnPfTYb0bURirbGoq9qmvCgznxriTzjFh3GKBo2tn5TievSWLbiYSaN3vVrfrOOUajR9CtrGY71dxwODOgu3nIauWOjo9Djt26EFrcL3W3ehhqUpZy9mn0eNzqD3oamnFmy2YhGWlvg9bhwLpYGoL2hMTpNyMxI5Bg9Ou5xX9hv6F6tZQbQ11vXm6aHmZHIs64RndH0K4oR1TqyqnQcmnxhXhgUlVR89W1arAUa7ypod+OqlNXIKa+EzUaUSNpbvliWFsa1LBIDQHebF0T6tzpIZAsosr7fdW/Yj/lMHsWSvuQZsVQOkdYW+Fu0avaBSH/jupQtoMT67nE07MdcWr8ThJTUUrvmnpCsuQnitLZCGBSVJPIMr9uFdhuzfZqV3VQvegLejGJ0GsUoFSNqk9swsPb96vmePW4Xutt8ld6+VvQujEvX+MAMzGf0LRTr9SBscbvQFfRhTueoTO/COCAZ0TIDCxmd37dOD0JJs1f3fbYL1X6CRPRWACPV1zDz31tQp6ZkKV9Gb9hn65a80bAPJyYTtpW3kUS+DH+LC2G/endSo0RDfpyaTtlW3kYS+TJaW9wIqcynZQaKo4feUVFvyIBBqeRq07EoH1pby9CTpsbIrpi9IZ/ueBC9C+NKuYA0XadnPTWW1h8c3RPy6zaidqFqhCK7CP9XAHcDuF1+jFpYr6YjkWPbPLwUekN+zGdWUdA5pWCUpZzkWmmrEY0Ym0YxiqTZ7o6DPELROZ1qZCRrpLduNAtuXE5lpIdo2Kd73aiy3bGuNRRjAZ1GttNu9OhdDWq7YaMA9nOjc2I0kGSeMTJg3zQIsObNMpfOY6Dd/nCfRJ5tnfoB1k+j2JGccSOJPKO33d5ylQVavR6E0kg2qetavYvEUrn6G9dymeVevr5OWjTsx0mdI9lK+nhdayj6jWi5zJhL53WPRKMhP043cPSuBrVrKCcB9FlZkWZnSUcglFHM3ExIDw0xKA1O1Z3Is20uwwrK70rvCLg35MfCcl7XSDaeziPYAs2LxADQFfTCRfrihpayqyiWWfVOjRsx4hCgaPZ5tGtWHAL0/D4XZc1GRmWNHL2rQa1B6QZwmogeJ6JHlYeVFWsmsqtFrBSt36lxI43cB4GZpekfm6f5lFFJI9wjK5ptv8/KGor+KS9pVKe91xxL5dDu0ze9pzgE6OmtG91nx4hDQCyVQ4dOzYpDgJ4MFkbWboAqhwADGQKsRu2U1x9YWYlmp7KXvI0uw0B1inD7501TuSJWy/Z6eAFrRrsRrsOplSIKDdB8284OfOT2Ybz1hi5d11dPw/RHtE2NxlJ53QYFUDL/6mhcDW6nXT2S1To1GkvnETGgWe9ahtKO9OhsR3qrgp3t/o2qRZVBYeaniSgKaTEeAF5i5rh11WoujMy5GqEz4IXHRQ0ZoVSimG2e8uoK+uB2UUOMqBEXWiMEfR586QO36L7eyFpGPJXDrqD+6IFoyI8ZHbERcQPOANXX6fFui6dyuKHNgGadThDKNXrXBiuam3hhXq2X14cBvAQpn9eHAbxIRB+0smLNhPKjtbuhcbmkoLdGTP9UpiRsnvJyuwg9bY3S3BiDYpSKK6vG70xZGO/wGxmh+HVO/yi9db3TP/rWF0uyZkOjMp1bBlQ068zF1SzJU7dCrZn+XQC3M/NHmfl/B3AHgM8ZLZyI7ieis0R0gYg+U+N9HxF9U37/RSIaqXrvs/Lxs0R0n9G6bEWlobF5ygtQEvDZ3yOZbWDjGm3Q7nRKFLLdayhG6WrzSYvjGnvryiKxkcZVWijW7toeS+XQGfTqWhgH9DsELCxLUe5GjejCsvbF8Vg6h+42L7wefaOjSlYEp49QALg2THEtaLi2JkTkBvAXAN4DYD+Anyei/RtO+ziAJWbeA+DLAP5UvnY/gI8AOADgfgB/KX+eJcTTebS4gHCrfcFuCn0NipZvZG89aiBozQiNGokaxe0i9IS055hSDKiRxlX5ruY0GjOjm9XpdQhQfldGjageh4BYUn/cDVCVFaGJ83mpNQo/kD28PkZEHwPw7wAeM1j2HQAuMPNFZl4F8A0AD2w45wEAX5effxvAO0mKOHsAwDeYOc/MlwBckD/PEhRPGDuD3RSi4cZM/8RTOQQ80JxbygyiBpMO6iWWyul2oW00evK+KVNVRqd/AO1TT0ai5BWiYb/m6R+lnnq9vAD9ru2xtPF9dnp1dBxOTCbwO98+gYUV692N1S7K/zYRfQDAT8uHHmbm7xksexDA1arXkwDu3OwcZi4SURJAl3z8hQ3XDtYqhIgeAvAQAESjUYyPj2uu6MxsDu3esq5rjbI8v4p0rojHf3QYPo99Bu3kxRzCXm6I5uzCKpIrBfzwycPwum3U/EYO4ZbGaDaKezWHi0ltdX/mqrTVsre0olvz1WQJAPDU80eRvKh+BH91PosOchv6rt2rObwxrU3zs7LmFgOaJ2XNT/7kCJai2jR3UtaQZk8hhzemM5o+4+nJAr55chUHRq3/bav+Npj5OwC+Y2FdLIGZHwbwMACMjo7y2NiY5s8YGwMOHz4MPdcaZSE0iW+fP459h+7ASHfQtnK/curH6GpNNUTzXNtVfOf8Cdx06E7s6ArYVu6XT/0YnaXGaDbKE0uv4fsnZzXV/dUfnQNOnUd/R1C35ng6hz94/kn07tiDsZ8aUXVNsVRG6vHv49CNIxgb26erXAD44dJr+OEpbZpfeeIc6LRBzSlJc8/OvRi7a6eqayqa9+3C2NiNusoFgMcXT+CJ03FNdT/+o/PAyXPo79SvWS319pR/Tv6bJqJU1SNNREZzAEwBGK56PSQfq3kOEXkg7aazoPJaU2nEdBdgPHeQXuKpHNp9jUlG3Shvlngqhw6/MxNwR8N+LC6vIl8sqb4mlsqju01yTdeLHjfvheVVlNn4ZnVRHbnu4ukcuoI+Y5oVJwgN/5PzmVWwzpT51ejJihBLSw4QLQY0q2XL/x5mvlv+G2LmcNUjxMxhg2W/DGAvEe0iIi+kRfaN0fePAvio/PyDAJ6S84k9CuAjshfYLgB7Ibk1bzsqmwnZaFDMcCc1gmJQ7Nz7oWyCO2kjUdYytCyOG0nOqKC4eWvp8Jjl8KG4DmvRLK3dGGvU3S6SHQJ0aDb4fevJihBL2hcIqSXbcN1jWmDmIoBPAHgc0h7132LmU0T0BSJ6v3za3wDoIqILAD4F4DPytacAfAvAaQA/APBrzKy+a+YgehsQzLSwbNyd1Ah9DRiVzcvupE41KHoC/cxYJJbK9iGmsVFXrjNarvR56n8nsyY1rtGwX9t3bZIR1ZOc0qz7rAa1aygHql/I00+3GS2cmR/DBm8xZv581fMcpGDKWtf+MYA/NlqHZifs96C1xW1r46qU1ajGNdzqgc/jsnUzIcVgN2pUZhQ9KWtiqTwODkQAZA2V3RPyY3JJ/WeYNkKpeFtpGJWlc3jzcATAsqGyo2EfphIavuu0OUZUT36/2aQ591kN9dZQPktEaQC3VK+fAIgBeMTy2glARLa7DivupI1qXCXNflunvMxwJ20ka2tt6hrXYqmM+Yxx112pbG0bfMVTObhICk40Vq4SA6Pud1IolTGfWTVFc6/GrarjqRzcLkKXzih5BcUgqf2+C6UyFpb1bQamh3prKH/CzCEAf7Zh/aSLmT9rSw0FtkfLx0wI/jKK3QGdFc0OHaFozfu2tkhszvSPFocAyRnAB4/bmANEV9CrySFgzsTA1WjIj4XlVawW1S2Ox1I59LRJDgxG0OoQMJfOgxm2bclQb4Ryk/z0n4noLRsfNtRPAH0BXEaYTeZABEMZWY3Sq7HXa5TZlKQ57HWmQVHyvqltXI2mUq8mqnFxXJrTN97ArWlW979hiWaVi+NmOAMA2rMimKlZDfXWUD4FKSjwv9V4jwG8w/QaCa6hL+zDE6kcmNkW92UzXCuNEg378dSZuH2aUzmp19xAzUbp0ZCosXodYz5mrNzequm2oY76cUOxVB6DJu2K2RtS7xCgGNvekB/zswbLrXIIGFSxm2oslcNwpzkxVVqyIph5n9WwpUFh5ofkv/daXxXBZkTDfuQKZaRWiogEWiwvb6031TjHub6wH9nVEtL5IsJ+OzTnGq7ZKNGQDxML6hablUa4N+zDvMFytWY7jqdyuHVHu8FS5bLDflxdVLfYHK/ansC4ZsX7UqXmdB6jIx0GS1XK9mFyaUXVudUbmRnVrAa1bsMfIqKQ/Pz3iOi7RHSrtVUTKFR6gDZNe5nlWmkEuzfamk3lG5JN2kw09VyT8iJx0IzpH/WeR6vFMhaWV037rqXNrtT9RmYrmo05A0jlqneCyBdLWDRRc68Gl+XZVA4ekzSrQe2q2OeYOU1EdwN4F6T4kL+yrlqCauyOy4jb6Le+GVq9lowST+Vs84Sxir6IH8mVArKrxbrnziSlbL9GF4kBySHA63ZhRsXvc22TKXN+X/2RVixlC8gV6o8sZ5M5REM+uEzQ3BX0osVNqjYXiyXlUYJJm9X1y04QqjWH/aZoVoNag6LU/L2QEkP+OwB7TJ5gLVreBjdaM10rjWBntLzSa7bLE8YqlLn8aRXxEdOJFVVz/2pwuQj97X5V5U4lpKmawXZz1hMG5LWY6UT9KaCpxAoGO0zUHGlVXS4ADJn0fQ/In6PGmE2ZeJ/VoNagTBHRVwH8HIDHiMin4VqBQSrztTZ4PZnpWmmESkSwDdN8iqdOo0dlRhmoGJT6jdx0cgX9JjY0AyobV+WcAZMW5fsjGoxocqXyHZnBQLtfk2azvm9N9zmxYtp3rQa1RuHDkFKk3MfMCQCdAH7bsloJ1tHqdSPs99gy5TVrs5vhZgS8HoT8Hlvib9Z2anT2CEVtb71cZswkcqY2NAPtWg2KOY3roMrGtVxmzCZzJhsUbZr7TZryUjRP1Sm7ZIHmeqgyKMycBfAGgPuI6BMAepn5h5bWTLCOPpu2xY2nmqdxjdoU3NhMmo0QDftBBEzXmQpZWF7Faqls6lTIQLt0r+ptizudlDLfmrWJ2ZrmrRvX+UwehRJjwKRGHZBGZbF0/a2Ap5PS1r+maZbXn2bqjMrmM3kUy2zqSLQear28PgngnwD0yo9/JKJft7JigvVIuxha31uvdjNsNH027dxod/CXVbS4XYiG6k/DVEYJEXN762VG3ZgQs6dgvB4XekO+upqnTB4ZKZ9VkrNUb4Wk2bxyfR43ejRoNivmRw1qp7w+DuBOZv68nLzxLgD/p3XVEmykN6Qtd5BeFDfDzkDjfS56w/bsLT+byqPFTehoAs1GUTOvvzanb+6UV/Vnb1W2mYYMgLw4vvX/hvK+2Wso0mfX12zWdNda2a11R2VmTy+qQa1BIayP+CrJxwQ20ReRUpGUy2xpObGU5E5ql5vhVkTlyG+rNSv7gjSDZqOomddf67ma19AMqmhcmRlTS+b21qWy62u2onFVs5bBzKaPUKSy/XXXUJrZoHwNwItE9AdE9AeQ9nP/G8tqJbiGaNiPUpkxv2xtjz2eypvmL2+UaMiHQomxmF21tBw794uwmsH2VkwntzbC04kcAl43Iq3mZSBQvK22auRSuSKWV0umu7EOyI2rtPdebaYSK2jzeRD2q98Dvh797fU9zFIrFmmWveq20jydyCHk89iSaUJB7aL8nwP4JQCL8uOXmPkrVlZMsJ61VA/WGpRYKtc0EeN9EXsCOqVUM82h2SgD7a2VuJrNUHrMZuZIC/o8aA+0bDlSsKrHPNDeinyxjMW6mv2mam7zeRBp3VqzFWs3yuflCmUsZQtblm3n6ASon23YT0S/QUT/A8DtAP6Smf87M79iT/UECnp2p9PDbKp5eut27VZp5xapVqM0IFtteGVVQzMQad0yx9TUkrkxKJVyK5q3btgt0dzeWve7Vs4zu1ygzn1esjcGBag/Qvk6gFEArwF4D4D/anmNBDVZ661b17hmV4tI54rNM+VlQ8qZ5XwR6Xxx2xiUkS4pAv3yQu2GhpkxMb9cOc/UsrsDm5YLoJK4cqQraG658udtlhhzTbO55Upl19E8v1w5z9Ryu6XPm9jqPi8sY6cFmreinkHZz8z/iZm/CuCDAO6xoU6CGnS3+UAES91oq9N7NwNKFltrNefWleV0hjsDIAIuzdduXBeWV5HOFy1qXIO4upjdNC7j0vwyIq0t6DA5UeFOubGemK/duM5l8lheLVlkRIO4spXmhWW0B1rQbrIH4c5O2Yhucp/n0nlkV0vY1d1cBqUyQcfM9TPOCSyjxe1CV9BnqeuwEjFutoujXlrcLnS3eS0dlSnGqlk0G8Xf4sZApBWXN+mtKw2QFQ3NSHcQxTJvujB/eSGLEQvKlTT7Nx2hKCMIK8re1SVp3mxh3qqRUavXjf4tNE9YqHkr6hmUN2/YS17ZWz5NRCk7KihYoy+iPlW3HmJNGDFudfxNRfM2MSiANB1yaZOpEGXkYkVDozScm42OLlk01QZIerYqt7p+ZqKMji5tYcCtGiXs7ApsOkKxaqqtHvX2lHdv2EveU/U8bFclBRLRkLXR8kr20r4malz7ItZGy1c0N5ERNcpIV3DzhmZhGW4XYcikrLvrylXm9WuUnSuUMJ1csaRRl8oObt5bn1+GxyLNirHYXHPOMs27uoObrqFcWpA025lpGGhQxmAi6iSiJ4jovPy35lZmRPRR+ZzzRPTRquPjRHSWiF6VH7321b5x9Iat762HfB60+czz1TeKtIGSdUY0lswh5Pcg2ESajbKrO4jkSgFLNdxoJxayGO5oRYvb/H/9njYfgl53zUbu6mIWzNZMtQHS1FMiW0CiRszSxMIyhjsD8FihOSRprjU6urKoTDtZNCrrCmJxeRXJlWtdhyfml7HDIs1b0agU9J8B8CQz7wXwpPx6HUTUCeD3AdwJ4A4Av7/B8PwCMx+SH3E7Kt1oomGflNivuHUyOr3MJFeaanQCSFNeC8t5FOok4NPLTDK3bdZPFLbyepqYt87zh4g2HSlYPadfWZivYcwm5rOWTf0QEXZ2BWuuWVk51QasfZe1yp6waL2qHo0yKA9AckmG/PfBGufcB+AJZl5k5iUATwC436b6NSXKtIyyf4fZzKbyTWdQomE/mKXMqVYg7SXfXJqNovSIN/aaFfdZKz1/Rrpqr2VYPaevaLo0n1l3nJlx2WL32ZHuwNaaLfq+N1uzWtNs7/oJADRqnB9l5hn5+SyAaI1zBgFcrXo9KR9T+BoRlQB8B8Af8SY5CIjoIQAPAUA0GsX4+LiuCmcyGd3XmkVsTnK0+/7hn2BPhzmpsKu5Es/iYLe7orMZNMfjkubHDv8EN7Sbr/nyXBZvajLNRimWGW4CfvTyaXSmLlSOz6+UsbxaQjkxjfHxucpxMzV7squ4slDA4z86DJ9nLSp9/EQeYS/h1Zd+Yko5G6lofuk0OpJrmueykma2UHNLdhWXFwp4/MnD8LnXND9zIo92H+GVF39sSjkbKcian3jpFCKJ85Xj8WwZWYs1b4ZlBoWIfgSgr8Zbv1v9gpmZiLRm//sFZp4iohAkg/KLAP6+1onM/DCAhwFgdHSUx8bGNBYlMT4+Dr3XmkXPdBJfPvocBvfsx9ib+k397GKpjOTj38ehfSMYG9sHoDk0d08l8ZVjz2FwzwGMHaz1c9JPQdZ8a5NpNoMbX3sWmRYfxsbuqBx74nQMwBE88PZR3LZzbfbYTM257lk88sZR9N54CLfuWCvjvxx/Fm/e6cXY2J2mlFOLvSeewbLXv07z46dmARzFA2OjeMsOqzTP4JE3jiF64604NNxeOf6nx5/Fm3euvwdms+f4M1huWa/5BydlzW8fXXcP7PhtWzblxczvYuaDNR6PAIgRUT8AyH9rrYFMARiuej0kHwMzK3/TAP4XpDWWbU+fhZHj85lVlLm5PLyANRfmuAVbAc+l82AG+kxOp94M3Nwfwusz6z37ldc39YUsK3d/f1guK105ViiVcSGeqbxnZdm1NBMB+6LWab65onmt7NViGRfi6cp71pUdWvddK/UgAvZZeJ83o1FrKI8CULy2PgrgkRrnPA7gZ4ioQ16M/xkAjxORh4i6AYCIWgC8D8BJG+rccDoCXrS4yRLX4Rl5b4Vmc5/tCnrhcVEl6NJMFHfkvsj2iJKv5uBABHPp/Lrv7cRkEru6g5Z6tA11tCLs9+C1qUTl2NnZNFZLZewfsNigDIQRS+XXeUK+NpnEri5rNQ93BBDye3BiMlk5dnY2jUKJLdd8cDCC2VRuveapJHZ3BxHw2r+i0SiD8iUA7yai8wDeJb8GEY0S0V8DADMvAvgigJflxxfkYz5IhuUEgFchjVr+p/0S7MflIssC/WKp5otBARTN1rgOz1ZiULbfCEWZ0jp6eQmAtFB77MrSuqkuK3C5CG/Z2VEpt7oOoyOdlpataDtSpfmoXZp3dOBYleYjlxcBAKMWl/0W+fOPXZHKLpel+zy609rvejMaYlCYeYGZ38nMe+WpsUX5+BFm/uWq8/6WmffIj6/Jx5aZ+TZmvoWZDzDzJ5m5tFlZ243esA8xC6Z/mjnAr1feaMtsZpswkNMs9g+E4W9x4eUJqWF7Y24Zi8urljdwgNSInotlKnEwL08soj/itzzI7sBABD6PCy9dUjRnkMgWMDpivebbdnbgXDxdiYM5MrGEgYjf8vTxByuaJYOiaL7NBs21aNQIRaCTvrDfmt56Kgev24VOkxP3mYEU3GjNlJfX40JHwL4NiOyixe3CXbu78OSZGJgZT74eAwD89J5uy8u+e28PAODw2ThWi2U8fW4Od9tQrtcjaX7qTBzMjB+9Li3N2qO5G8xrmp85N4e799qj+c7dXXhKvs+KZju+71oIg+IwomE/YhasJ8SSOUQjPlM3IDKLvrDfmjWUZA59YXM3XWom7j/Qh6uLK3j1agL/dmIGBwbCGO60PjbhlsEI+sJ+/OvxaTxzbg7pXBH3HTDXQ28z7j/YhyuLWZyYTOJfj0/jTYMRDHVYr/nQUDuiYR/+9fgMnj43h3S+iPtN9krcjPsORDGxkMVrU5LmW4Yitm+spSAMisPoDfuQzhexnDc3+fOM3Lg2I71hP1K5IlZWzZ3ZnG1izWbws7f0I9Lagv/wlz/Ba1NJ/MKdO20p1+Ui/Mc7d+Dw2Tk89A9HMNTRintu7LGl7PfKmh/4ix/j1HQKv3DnDlvKdbkI//GOnXjqTBwP/cMRDHe24u499mh+3y0DCPs9eP//+DFOz9inuRbCoDiMvoobrbnTXrFUrmndZ61yHZ5N5bbl+olC2N+CLz54EP4WF955Uy8+NDpkW9kfv3sX7hjpRNDnwR89eBBejz1NTdjfgi88cAD+FhfedXMvPnCbjYJeI1QAABGnSURBVJrftgu3j3SgzefBHz/4Jts0R1rX7vO7bo7iA2+xT/NGtk9GvOsEpXGdTeZMS6HBzJhJ5vDu/c3pPttXpdmsFBrMvO0NCgC8/80DuO9AFD6P+VkGtiLo8+Cbv3IXVktl28t+4NAg7j/YZ3u5bT4PvvUrP4VCiW0zJgqN0rwRYVAchrLfu5m99eRKAfliuYlHKJLmmImjsqVsAavF8rae8lJoVCNDRA0ru5GavZ7GrMk12pgAYsrLcfRaEC3fzC7DwJpmM+NvtrPLsEDQKIRBcRghnwcBr9tU1+HZJg1qVAj7PWhtcZtqRGdTcmaAJtUsEDgRYVAcBhEhGjZ3F8Nm761Lmn2mppyZTUqf1ayjMoHAiQiD4kB6Qz7Tp3+IpM9tVnrDfnNHKMkVuEjacU8gEJiDMCgOJGpytPxsMofuNp8l28KaRdTk7Y9nHKBZIHAa4r/JgfRFpCmvTfYU08xUYsXyPEtG6Qv7zNfc0dyaBQKnIQyKA+kN+bBaLCO5UjDl85zQuEbDfuQKZaRy5mQIcIIRFQichjAoDiRacR02Pu1VLjOmEisYavLG1UzX4XKZMe0AIyoQOA1hUByI4o1lxiL1fCaP1WIZQ03euPaZaETj6TwKJbYlaaBAcD0hDIoDiYbWUpEYZTIhxWM0e29diZZXdpY0wlQiCwBNPyoTCJyGMCgOpC/iB9GaMTDC1JJsUNqbu7euaJ4yQfPkkjOMqEDgNIRBcSBejwt9YT8mF7OGP2vKISMUn8eNaMiPq4tmjFAUI9rcmgUCpyEMikMZ7ghUetpGmFpaQaS1BW2+5s8TOtzZisklE4zo0graAy0IOkCzQOAkhEFxKEOdrbhqRuPqIPdZ04yogzQLBE5CGBSHMtQRwGwqh9Vi2dDnTC05x312qKMVM8kVFEomaBYGRSAwHWFQHMpwRyuYgWkDi9TM7Kje+lBnAGUGZhL6vdsqmh1iRAUCJ9EQg0JEnUT0BBGdl/92bHLeD4goQUT/tuH4LiJ6kYguENE3ichrT82bh+FOySvLyLRXaqWITL7Y9DEoCsMdxjUnsgVkV0siBkUgsIBGjVA+A+BJZt4L4En5dS3+DMAv1jj+pwC+zMx7ACwB+LgltWxiFCNgZE1BaZidYlCUel414N1WcRl2yKhMIHASjTIoDwD4uvz86wAerHUSMz8JIF19jIgIwDsAfLve9duZ/kgrPC4y1Liuuc86o7feH/HD7SJDRrQS1OgQIyoQOIlG+U1GmXlGfj4LIKrh2i4ACWZWsgROAhjc7GQiegjAQwAQjUYxPj6uvbYAMpmM7mutosMHHDs7gXH/rK7rxyek5JITp49h4cK1+2A3o+ZOH3DkzCWM+2bqn1yD8Uuy5lPHMH/eGZqtRmi+PrBDs2UGhYh+BKCvxlu/W/2CmZmIzMlJXgNmfhjAwwAwOjrKY2Njuj5nfHwceq+1ir3nX8BKoYSxsZ/Wdf3hR04i5JvC+949Bmngt55m1Lzn3AvIF/VrfjJxEiH/FN7rIM1WIzRfH9ih2TKDwszv2uw9IooRUT8zzxBRP4C4ho9eANBORB55lDIEYMpgdR3JcEcAT57R8tWt5/JiFju7AzUb1mZluLMV42fndF9/eTGLka6gozQLBE6hUWsojwL4qPz8owAeUXshSzssHQbwQT3XbyeGO1sxn8kjVyjpuv7yQhY7O4Mm18pahjsCiKf1a76ysIydXc5YMxIInEajDMqXALybiM4DeJf8GkQ0SkR/rZxERM8C+GcA7ySiSSK6T37rdwB8ioguQFpT+Rtba98kKK6vetKRFEtlTC5lHde4DnXq926TNK84TrNA4BQasijPzAsA3lnj+BEAv1z1+m2bXH8RwB2WVdAh7JAbxssLWezpDWm6diaZQ6HEjmtcd8gjqiuLy9jT26bp2ulEDsUyY2eXs0ZlAoFTEJHyDmZ3t9QwXpxb1nztxIJ0jdMa1xt6TNDc6SwjKhA4BWFQHEx7wIvOoBcX57U3rpcXpGkyp41Q2gNedARa9GmWY3ZGup1lRAUCpyAMisPZ1R3EpfmM5usuLyzD63FVdn90Eru6g7ikY4RyeX4Z/hYXekM+C2olEAiEQXE4kkHR3rheml/Gzs4AXC7nuc/u6m4zoFm4DAsEViEMisPZ1R1ELJXHcr5Y/+QqLsQzmhe1m4XdPUHMpnLaNc85V7NA4ASEQXE4ysK8lh57rlDClcUs9jq0cd0la1YW2dWQK5RwdTErDIpAYCHCoDicXT3aDcrEwjLKDNzg0MZ1lw4jemle0iwMikBgHcKgOBwpjYg0haUW5VynNq7Xo2aBwAkIg+Jw/C1ujHQFcXY2Xf9kmQvxDIiAG3qc2bi2et3Y2RnQrNlFa6MbgUBgPsKgbAP2RUM4F9PWuA51tMLf4rawVtayry+Es1o0z2Uw3BlwtGaBoNkRBmUbsK8vhImFZdUJEy/EM9jj0NGJwr6+MCbm1Wt+YxtoFgiaHWFQtgH7+kIoM3A+Vn9NIV8s4UI8g319YRtqZh37opJmNesoa5q15TsTCATaEAZlG6A0lGdmU3XPPR/LoFhmHBx0uEGpaK4/7XVuVtEcsbpaAsF1jTAo24CRriC8HpeqdZRT00kAwIEBZzeuI10BHZqdbUQFgmZHGJRtgNtFuKkvhJNT9Ucop6ZTaPN5HJ9x1+N2YV80hJNTybrnnpxOIuTzYLjD2ZoFgmZHGJRtwqHhdpyYTKBU5i3POzmVxP7+sCNzeG1E0pysq/nUdAr7B7aHZoGgmREGZZtwaLgdy6ulLRepS2XG6zNp7N8mUz+HhtuRyRfxxlw9zSnHT/EJBE5AGJRtwq07OgAAr1xZ2vScs7NprBRKuGVoezSut+5oB7C15tdnUsgVyttGs0DQzAiDsk0Y6QqgPdCCV68mNj3nxUsLAIA7d3fZVS1L2dUdRKS1nuZFAMCduzvtqpZAcN0iDMo2gYhwaLgdx7borb9wcQFDHa0YbG+1sWbWUdF8eQuDcnEBOzoD6I9sD80CQTMjDMo24q7dXTgXy2A2mbvmvXKZ8dKlRdy1TUYnCnft7sLZWHpzzROLuEuMTgQCWxAGZRvxjpt6AQDjZ+PXvHcunsZStoA7d22vxlXR/PS5azWfjaWRyBZw567tZUQFgmalIQaFiDqJ6AkiOi//7djkvB8QUYKI/m3D8b8joktE9Kr8OGRPzZubvb1tGGxvxeEaBmX87BwA4K17uu2ulqXcGG3DQMSPw2fmrnlP+R7eukcYFIHADho1QvkMgCeZeS+AJ+XXtfgzAL+4yXu/zcyH5MerVlTSaRARxvb14Lnz88gX1ydN/MHJWbxpMLJt1k8UiAhjN/XiuQvXan785CzePBQR6ycCgU00yqA8AODr8vOvA3iw1knM/CQA9TnKBbjvQB+WV0t44nSscuxCPINXrybw3lv6G1gz67jvQB8y+SJ+dHptZHY+lsbxyeS21SwQNCPEvHWUsSWFEiWYuV1+TgCWlNc1zh0D8FvM/L6qY38H4KcA5CGPcJg5v8n1DwF4CACi0eht3/jGN3TVOZPJoK2t+dOfl5nx20+voNNP+M93+kFE+PvTeTx9tYg/Hwsg4lMfLe40zV2thM/eIWn++qk8npks4stjAYS3oWYzEZqvD4xovvfee48y82i98zy6Pl0FRPQjAH013vrd6hfMzESk1ap9FsAsAC+AhwH8DoAv1DqRmR+Wz8Ho6CiPjY1pLEpifHwceq+1m99ovYzP/ctJpDpuxN7eEJ594jl8+PYdeOC+N2n6HCdp/qR/Ap9/5BQynfuwuyeI5574MT5yxw68fxtrNguh+frADs2WGRRmftdm7xFRjIj6mXmGiPoBXLuKvPVnz8hP80T0NQC/ZaCq246P3D6M7x2bxG9+81V43C50Br341LtvbHS1LOUjt+/A916Zwie/8Qo8Lhe6gj785jbXLBA0G5YZlDo8CuCjAL4k/31Ey8VVxoggrb+cNL+KzqXF7cLXPnYH/t+nziOdK+LX7t2DnpCv0dWyFK/Hhb/72B3470+dR0bW3N22vTULBM1GowzKlwB8i4g+DuAygA8DABGNAvhVZv5l+fWzAG4C0EZEkwA+zsyPA/gnIuoBQABeBfCrDdDQ1EQCLfi99+1vdDVsJRJoweeuM80CQTPREIPCzAsA3lnj+BEAv1z1+m2bXP8O62onEAgEAj2ISHmBQCAQmIIwKAKBQCAwBWFQBAKBQGAKwqAIBAKBwBSEQREIBAKBKQiDIhAIBAJTEAZFIBAIBKbQkOSQjYKI5iAFUuqhG8C8idVxAkLz9YHQfH1gRPNOZu6pd9J1ZVCMQERH1GTb3E4IzdcHQvP1gR2axZSXQCAQCExBGBSBQCAQmIIwKOp5uNEVaABC8/WB0Hx9YLlmsYYiEAgEAlMQIxSBQCAQmIIwKCogovuJ6CwRXSCizzS6PlZDRH9LRHEium42LiOiYSI6TESniegUEX2y0XWyGiLyE9FLRHRc1vyHja6THRCRm4heIaJ/a3Rd7IKIJojoNSJ6lYiOWFaOmPLaGiJyAzgH4N0AJgG8DODnmfl0QytmIUR0D4AMgL9n5oONro8dyFtR9zPzMSIKATgK4MFtfp8JQJCZM0TUAuA5AJ9k5hcaXDVLIaJPARgFEGbm9zW6PnZARBMARpnZ0tgbMUKpzx0ALjDzRWZeBfANAA80uE6WwszPAFhsdD3shJlnmPmY/DwN4HUAg42tlbWwREZ+2SI/tnUPk4iGALwXwF83ui7bEWFQ6jMI4GrV60ls84bmeoeIRgDcCuDFxtbEeuTpn1cBxAE8wczbXfNXAHwaQLnRFbEZBvBDIjpKRA9ZVYgwKAJBFUTUBuA7AH6DmVONro/VMHOJmQ8BGAJwBxFt2ylOInofgDgzH210XRrA3cz8FgDvAfBr8rS26QiDUp8pAMNVr4fkY4JthryO8B0A/8TM3210feyEmRMADgO4v9F1sZD/v737CbGqjMM4/n1CpND+gAxSuDCssVykgoYwhBpR4KIWbYKxVZugPxS4sVWLFgoy0cqVu8QQ7B8MlUFRNAVNUDqVbmqlGA7RooFJCp8W5x08t8aZiPfcw+Dzgcs977kv5/zu4vK773nf8ztjwONlPuEt4GFJb/Yb0nDYvljeLwPv0FzKry4JZXnTwL2S7pa0GngKeL/nmKKyMkF9DDhne6LveIZB0oikO8r2LTQLT873G1V3bB+0vcH2Rprf8Se29/ccVuckrSkLTZC0BngU6GQFZxLKMmz/BTwPfEQzUXvS9g/9RtUtSSeAr4DNki5IeqbvmIZgDHia5l/rd+W1r++gOnYn8KmkszR/nD62fcMspb2BrAe+kHQG+BqYtP1hFyfKsuGIiKgiI5SIiKgiCSUiIqpIQomIiCqSUCIioooklIiIqCIJJSIiqkhCifgfJK1r3a/yi6SLrfaXHZ1zu6RjS3w+IqmT+wsi/otVfQcQsRLZ/hXYBiDpVWDO9pGOT/sK8NoSMc1KuiRpzPZUx7FE/EtGKBGVSZor73skfSbpPUk/Szokabw81GpG0qbSb0TSKUnT5TW2yDFvBR6wfaa0d7dGRN8ulNYA3gXGh/RVIwYkoUR0ayvwLHA/TWmXUdsP0jyP44XS5w3gdds7gSdZ/FkdOxisv3QAeK5UCn4ImC/7vyntiKHLJa+Ibk3bvgQg6SfgdNk/A+wt248AW5r6lADcJmlt6+FX0NTdmm21p4AJSceBt21fKPsvA3fV/xoRy0tCiejWldb21Vb7Ktd+fzcBu2z/scRx5oGbFxq2D0maBPYBU5Ies32+9Jm/zjEiOpVLXhH9O821y19I2rZIn3PAPa0+m2zP2D5MUyn4vvLRKB2VJo9YThJKRP9eBHZIOivpR5o5lwFl9HF7a/L9JUnfl9LzfwIflP17gclhBB3xTylfH7FCSHoZ+N32YpP2C30+B56w/dvwIotoZIQSsXIcZXBOZoCkEWAiyST6khFKRERUkRFKRERUkYQSERFVJKFEREQVSSgREVFFEkpERFTxN3cL1cvo6mtmAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "# simulation timing\n", "Tf = 5\n", "Ts = 0.01\n", "N = int(Tf/Ts)\n", "\n", "# initital condition\n", "x = 0\n", "\n", "# simulation history\n", "x_hist = np.zeros((N,1))\n", "\n", "for i in range(N-1):\n", " # dynamics\n", " freq = 1 if i < N//2 else 2\n", " v = np.cos(2*np.pi*freq*(Ts*i))\n", " f = lambda x: v\n", "\n", " # propagate dynamics by solving the differential equation\n", " x = rk4(f, x, Ts)\n", " \n", " # add to history\n", " x_hist[i+1] = x\n", " \n", "plt.plot(np.arange(0, Tf, Ts), x_hist)\n", "plt.grid(); plt.xlabel('Time (s)'); plt.ylabel('Position (m)')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Lie Group Integration\n", "\n", "Note that the above formulations we assume Euclidean spaces instead of matrix Lie groups." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.5.2" }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }