jax_for_hamiltonian.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Copy of jax_for_hamiltonian.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "PZkO10ct_0Ob"
      },
      "source": [
        "import jax.numpy as jnp\n",
        "import jax.ops as ops\n",
        "from jax import grad, jit, vmap\n",
        "from scipy.integrate import solve_ivp\n",
        "import matplotlib.pyplot as plt"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QSYmtqGI_5zU"
      },
      "source": [
        "# Learn `jax` basics"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "f3MWdWzQIhSN"
      },
      "source": [
        "## Differentiate a built-in function using `grad`:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "84vZbgqxWCpS",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "52696441-8a32-4127-973c-d081fd621757"
      },
      "source": [
        "# Differentiate cosine\n",
        "my_sin = grad(jnp.cos)\n",
        "x = jnp.pi/2\n",
        "y = jnp.cos(x)\n",
        "dydx = my_sin(x)\n",
        "\n",
        "print(f'y = cos(x)\\nfor x = {x:.2f}:')\n",
        "print(f'\\ty = {y:.2f}\\n\\tdydx = {dydx:.2f}')"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "y = cos(x)\n",
            "for x = 1.57:\n",
            "\ty = -0.00\n",
            "\tdydx = -1.00\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QRH1mL1zIrET"
      },
      "source": [
        "## Define my own function and differentiate it:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "VDCw9fx8IdOI",
        "outputId": "b783ba86-3d57-4441-95ad-20e79d091bdc"
      },
      "source": [
        "# Differentiate 2nd order polynomial, scalar\n",
        "def my_quad(x, p):\n",
        "  assert len(p) == 3\n",
        "  return p[0] + p[1]*x + p[2]*x**2\n",
        "\n",
        "x = 1.0\n",
        "p = jnp.array([1.0,2.0,3.0])\n",
        "y = my_quad(x,p)\n",
        "my_lin = grad(my_quad)\n",
        "dydx = my_lin(x,p)\n",
        "\n",
        "print(f'quadratic equation, y = p0 + p1*x + p2*x^2\\nfor x = {x:.2f}, p = {p}:')\n",
        "print(f'\\ty = {y:.2f}\\n\\tdydx = {dydx:.2f}')"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "quadratic equation, y = p0 + p1*x + p2*x^2\n",
            "for x = 1.00, p = [1. 2. 3.]:\n",
            "\ty = 6.00\n",
            "\tdydx = 8.00\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Jq4QOItmJF97"
      },
      "source": [
        "## Vectorize my function and its derivative using `vmap`"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "PvaOki6hId4d",
        "outputId": "8b0cd9c8-3cf6-497e-bce8-c8b58c760cf6"
      },
      "source": [
        "# Differentiate 2nd order polynomial, vector\n",
        "x_vec = jnp.arange(0.0, 1.2, 0.2)\n",
        "y_vec = vmap(my_quad, in_axes=(0,None))(x_vec, p)\n",
        "dydx_vec = vmap(my_lin, in_axes=(0,None))(x_vec, p)\n",
        "\n",
        "print(f'Using y and p as defined above\\nfor x = {x_vec}:')\n",
        "print(f'\\ty =    {y_vec}')\n",
        "print(f'\\tdydx = {dydx_vec}')"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Using y and p as defined above\n",
            "for x = [0.  0.2 0.4 0.6 0.8 1. ]:\n",
            "\ty =    [1.        1.52      2.28      3.2800002 4.52      6.       ]\n",
            "\tdydx = [2.        3.2       4.4       5.6000004 6.8       8.       ]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GxL3hwqFs9_q"
      },
      "source": [
        "# Solve Hamilton's equations over for a simple mass-spring system, using `jax`\n",
        "- $N$ masses of mass $m$ connected by springs with stiffness $k$\n",
        "- Right boundary condition: free ($N-1$ springs) or fixed ($N$ springs)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nu64uFtESIgR"
      },
      "source": [
        "## Define $V$, $T$, $H$ (using `vmap`)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "aUK1YD9q87An"
      },
      "source": [
        "def V(dx, k):\n",
        "  ''' Potential energy of deformation of spring between masses n and n+1,\n",
        "      with spring constant k and deformation dx = x_n+1 - x_n '''\n",
        "  return 0.5*k*dx**2\n",
        "\n",
        "\n",
        "def T(p, m):\n",
        "  ''' Kinetic energy of mass m with momentum p '''\n",
        "  return 0.5*p**2/m\n",
        "\n",
        "\n",
        "def H(x, p, params):\n",
        "  ''' Given positions x and momentums p, compute Hamiltonian of system with N masses.\n",
        "      Dictionary params includes stiffness k, mass m, and boundary condition BC\n",
        "      defines right edge of chain '''\n",
        "\n",
        "  if params['BC'] == 0: # free\n",
        "    dx = jnp.pad(jnp.diff(x),(0,1),'constant')\n",
        "  elif params['BC'] == 1: #fixed\n",
        "    dx = jnp.diff(jnp.pad(x,(0,1),'constant'))\n",
        "  else:\n",
        "    assert False, 'Specify right boundary condition as 0 (free) or 1 (fixed)'\n",
        "\n",
        "  V_vec = vmap(V, in_axes=(0,None))(dx, params['k'])\n",
        "  T_vec = vmap(T, in_axes=(0,None))(p, params['m'])\n",
        "  H_total = jnp.sum(jnp.add(V_vec, T_vec))\n",
        "\n",
        "  return H_total"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2-H2D4g7NX38"
      },
      "source": [
        "## Differentiate Hamiltonian and define equations of motion (using `grad`)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hIQIhc8wOYVl"
      },
      "source": [
        "dHdx = grad(H, argnums=0)\n",
        "dHdp = grad(H, argnums=1)\n",
        "\n",
        "def eom(t, y, params):\n",
        "  ''' Equations of motion for general 1D Hamiltonian with N units '''\n",
        "\n",
        "  x, p = jnp.split(y, 2, axis=0)\n",
        "  dxdt =  dHdp(x, p, params) # velocity\n",
        "  dpdt = -dHdx(x, p, params) # force\n",
        "  dydt = jnp.concatenate((dxdt,dpdt), axis=0)\n",
        "\n",
        "  return dydt"
      ],
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oWtXqSM93tlR"
      },
      "source": [
        "## Set chain length, parameters, and initial conditions"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QtsJmUTh3spr"
      },
      "source": [
        "def generate_IC_vector(x0, p0, N):\n",
        "  ''' Given inital position x0 and momentum p0 of the leftmost unit of an N-unit\n",
        "      chain, return y0 vector containing initial x and p vectors assuming the\n",
        "      remaining masses are at rest in equilibrium.\n",
        "      \n",
        "      For use in solve_ivp, need y0.shape = (2N,)'''\n",
        "\n",
        "  x = jnp.zeros(N)\n",
        "  p = jnp.zeros(N)\n",
        "  x = ops.index_update(x, 0, x0)\n",
        "  p = ops.index_update(p, 0, p0)\n",
        "  y0 = jnp.concatenate((x,p), axis=0)\n",
        "\n",
        "  return y0"
      ],
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Aj2wQrOhITPb"
      },
      "source": [
        "N = 25 # chain length (number of masses)\n",
        "params = {'k':1, 'm':2, 'BC':1} # {k: spring constant, m: mass, BC: [0|1]}\n",
        "x0 = 0.0 # Initially in equilibrium position\n",
        "p0 = 1.0 # Give initial \"kick\" to left end\n",
        "\n",
        "y0 = generate_IC_vector(x0, p0, N)"
      ],
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SVuPTlnn8oiG"
      },
      "source": [
        "## Integrate with `scipy.integrate.solve_ivp`"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Gne80J5R8vh2"
      },
      "source": [
        "# Set the solver parameters and simulation time\n",
        "tspan = (0.0,100.0)\n",
        "method = 'RK45' # Typical Runge-Kutta method\n",
        "rtol = 1e-8\n",
        "atol = 1e-8"
      ],
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "3cshICt0OY_i",
        "outputId": "a7b341c6-47a8-43cb-d5be-aa56e2065517"
      },
      "source": [
        "# Solve equations of motion \n",
        "sol = solve_ivp(eom, tspan, y0, method=method, rtol=rtol, atol=atol, args=(params,))\n",
        "\n",
        "print('Solution found over {0} time steps'.format(len(sol.t)))"
      ],
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Solution found over 687 time steps\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vKOC1vruOqJN"
      },
      "source": [
        "## Inspect results"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pasRUmG_RD2l",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 545
        },
        "outputId": "be2a92fe-ebcd-4822-bb6a-88a3449db33f"
      },
      "source": [
        "# Plot displacement of first mass over simulation time\n",
        "plt.figure('solution')\n",
        "plt.plot(sol.t, sol.y[0])\n",
        "plt.xlabel('Time (s)')\n",
        "plt.ylabel('Position of first mass')\n",
        "\n",
        "# Find and plot spring strain for all springs over time\n",
        "x_arr, p_arr = jnp.split(sol.y, 2, axis=0)\n",
        "dx_arr = jnp.diff(x_arr, axis=0)\n",
        "\n",
        "UNIT, TIME = jnp.meshgrid(jnp.arange(N-1),sol.t)\n",
        "plt.figure('colorplot')\n",
        "plt.pcolor(UNIT, TIME ,jnp.transpose(dx_arr),cmap=plt.cm.jet)\n",
        "plt.xlabel(\"Unit index\")\n",
        "plt.ylabel(\"Time (s)\")\n",
        "cbar = plt.colorbar()\n",
        "cbar.set_label('Strain')"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZhcdZno8e9bVV1dvW/pLd2dhSyEhJAAzSagyCaMCIjIhdERFR+u9+pl1HFGvHPHO+O93gdxZvTqoFdGURwdQHEhjiwii6AYSAeyk5Cmk5Du9Jb0nl5qe+8fdSpUOt3pquqqrqrO+3meeqrO75xT560U1Nu/5fx+oqoYY4wxiXJlOgBjjDG5yRKIMcaYpFgCMcYYkxRLIMYYY5JiCcQYY0xSPJkOYC4tWLBAlyxZkukwjDEmp2zevPmwqlZPLj+lEsiSJUtoaWnJdBjGGJNTROTAVOXWhGWMMSYplkCMMcYkxRKIMcaYpFgCMcYYkxRLIMYYY5JiCcQYY0xSLIEYY4xJiiWQFNjePsi//Wk/g2OBTIdijDFzJqM3EorINcD/BdzA91T1nkn7vw6829ksBGpUtdzZFwK2O/veUtXr5ybq47X2DPOB77yEPxTmZ5vb+cV/eQcet+VlY8z8l7FfOhFxA/cB1wKrgdtEZHXsMar6WVVdr6rrgW8Bv4jZPRbdl6nkAfCd59vwelz8j/eewbb2QR7edDBToRhjzJzK5J/K5wOtqtqmqn7gYeCGkxx/G/DQnEQWJ38wzJM7OvmztXXccclSzl9ayTd+9wbjgVCmQzPGmLTLZAJpAGL/XG93yk4gIouBpcCzMcU+EWkRkY0icuN0FxGRO53jWnp7e1MR9zHb2gc46g9x+aoaRITPXrmSwyN+fv5qe0qvY4wx2ShXGutvBR5V1dg/7RerajPw58A3RGTZVCeq6v2q2qyqzdXVJ0wmOSub9vcDcMHSKgAuPK2SdY1lfO/FfYTCJ19rPhgKs7triNaeYcIzHGuMMdkok53oHUBTzHajUzaVW4FPxRaoaofz3CYizwNnA2+mPszpvdE9zMIyHxVFXgBEhDvfuYxP/furPL2ri2vOrJ/yvJfbjvDXj27jrb5RAJYuKOLzV5/On62tQ0SmvV4gFGZP1zA9w+NUF/tYWVdMvsed+g9mjDFxyGQC2QSsEJGlRBLHrURqE8cRkVVABfCnmLIKYFRVJ0RkAXAxcO+cRB1jT9cwK2pLjiu75sw6FlUW8t0X2njPmhMTwiv7+vjoDzZRV+bjHz+4jlA4zA/+uJ9P/furvGtlNf/rhjNZVFV43DndQ+P88KX9/GTjAYbGg8fKfXku3rFsAe9eVcPlq2poKC84tm9wNMAr+/vY2HaEjW1H6BgYIxRSakrzWV5TzKq6Us6oL+H0ulIWVxYiAp2D47zRPUxrzwh7u0doHxhlPBDG63ZRW5pPQ0UBi6uKWFJVRENFAXkuYSwQ4tDAOIcGxiKPwTEGxwKEw1Ds81BTkk91ST41JT5qSvMp9eURCIWZCIYZHPMzMBqgfzTA4Kif/tEAwXCYsgIvFYV5VBR6KS/Mo7zQi9sl+INhxoMhRsaDjEwEGRkPMuw857mF0oI8ypxHoddNQV4kuY76Qxz1Bzk6EeLoROTcUX8Qj9tFYZ6bAq+bonwPhV43hV4PXo+LiUCI8WCYcX+IsYDz8Ifwh8L4PO7I+3vd+PLceD0uvG4XeW5hPBBm1B9k1B9i1B8iEAqjCooiCHkeIc/lwuMWPG4XeS7n2S0EQ4o/FCYQCuMPhvGHwoRCissleFyCO/rsdh3bdokQdM4JhJRAKEzQqdG6RHC7Is+R14JLwOUS3CKEVAmHlWBYCYUVlwhF+ZF/g2PPXg+lBR5KfXmU+Dw2wtAcJ2MJRFWDIvJp4Ckiw3gfUNWdIvJloEVVNziH3go8rKqx7TxnAN8VkTCRZrh7VHXXXMYfCiutvSNcsmLBceVul/CJS5fypcd20nKgn/OWVB7b99pb/Xz8h5tYWO7j4TsvorokH4APnNPIj/50gH/67R6u+vrv+a+XLeey06vpG/WzYcsh/mPbIUJh5T1r6rh2bT2LKgvp6B9j0/4+nt3dw7O7e/g7YGVtMVVF+XQNjbP/yFFUwetxce6iCs5dV4FLhM7BMfb2jPD0rm6iLWd5biGsHNfsVlXkZVFVIYVeNxOBMJv29/PrbZ0zNs0tKM6nvDAPl8DQWJDDIxPHftBORgTKCvLwuFwMjvkJhOJv1ivyugmEFX8wHPc5vjwXwZDGFVssj0sSPmc+KfS6KfXlUVOazxl1pZzVVMaVZ9RSW+rLdGgmA+T43+X5rbm5WVO1oFRb7wiX/9Pv+drNZ/HB5qbj9o35Q7zjnmc4rbqYR+68EI/bxY6OQf78XzdSXujlp//5IurKTvwfrmtwnC//x04e3951rKw438PN5zbysYuXsLiq6IRzVJU3e4/y3O4eXmw9zJg/SGWRlzMXlnH+0krWNZXjyzuxmWvMH2JvzzC7O4d58/AIHpdQV1bAyppiltcUU1Wcf8I5gVCYjv4xDvSNcmhgjLAqXreLheUF1Jf5WFhecMK1wmGlb9RPz9AEPcPjHJ0IkecWvB4X5YVeygsiNY0SnweXS459pqP+EP1HozUUf+RaHhf5HjelPg/FPg/F+ZG/kKPnjQdCDI4FGBwLMO7UGBQo8noo8Lojxzt/Wbudc/zBMGP+EKOBSA1lzB9iIhjCl+d2Hi4KvR58eS58HjcuV+Qv/miNZDwQxh8KMREMEwwpvjy3U5OJXCfPLYgIAqjzbxgIRY4NhCPP0dd5Lhd5HonUZpxajccVqSmEojWF0Ns1hmA4TDgMHreQ53aOdwsetzj/9hBWPVbTiP6RENbII1or8bgEl0sIhZTRQIjRiSBH/ZHnkYkgw+NBhsYDkeexAEPjAQ4NjLPz0CD9o5GbZy88rZJPvmsZ71pZfdJmWJObRGSz0+d8fLklkOQ8uaOLT/54M4996mLWNZWfsP9Xr3XwmUe2cNM5DZzdVM49T+ymrCCPn37yIhorCqd4x7ftP3yUN7qHKfHlsa6pjELvKbVwpMkRkT9eRnh8excPvfIWnYPjnLekgv/z/rUnNO2a3DZdArFfpiS190c6wJdMUSsAuPHsBlp7Rrjv+VZ+8WoHZy8q5zsfOnfKmsdkSxYUsWTB1O9rTLYQEZbXlHDXFSV88l3L+Nnmg9z75B7e+60/8H/ev5abz23MdIgmzSyBJKljYIwir5vSgun/CT//ntP50IWLGBwLcHptiVXtzbzl9bj40AWLuXp1HXc99Bqf/9lWDg2McdcVKzIdmkkjG1KRpEMDYywsL5gxKdSXFbCqrtSShzklVJfk8293nM9N5zTwz0+/wbefb810SCaNrAaSpEMD4yyMGTZrjInwuF187eZ1BELKvU/uYWVNCVeurs10WCYNrAaSpGgNxBhzIrdL+NrNZ3FmQymf/ekWuofGMx2SSQNLIEkYD4Q4ctRPQ7mNfTdmOr48N/9y2zn4g2G+9NiOTIdj0sASSBIODYwBWA3EmBksWVDEZ69ayVM7u3l2d3emwzEpZgkkCV2Dkep4PENyjTnV3XHJUhZXFfKPT71hE4fOM5ZAktA7MgFATcmJd2sbY46X53bxmStXsKtziCd2dM18gskZlkCScHjED0TmfTLGzOz6dQ2cVl3Ed1+Y0wmzTZpZAknC4ZEJ8txCWUFepkMxJie4XcJH37GEbe2DbDk4kOlwTIpYAknC4eEJqory7eZAYxLw/rMbKPK6+dGf9mc6FJMilkCS0DsywYISb6bDMCanlPjyeP85DfzHtk6GxwOZDsekgCWQJBwemaDa+j+MSdj7z27EHwzz9C4b0jsfWAJJwuFhv3WgG5OEcxaV01BewK+3Hsp0KCYFLIEkSFU5cnSCBTaE15iEiQjXravnxb2H6T/qz3Q4ZpYsgSRocCxAIKRWAzEmSdetXUgwrDyzuyfToZhZymgCEZFrRGSPiLSKyN1T7P+oiPSKyBbn8YmYfbeLyF7ncftcxdzn/NVUWWRDeI1JxpkNpVSX5PP8HksguS5j07mLiBu4D7gKaAc2icgGVd016dBHVPXTk86tBP4n0ExkqenNzrn96Y57YCwyeqS80EZhGZMMEeGyldX8dlc3wVAYj9saQnJVJr+584FWVW1TVT/wMHBDnOe+B3haVfucpPE0cE2a4jzO4GgkgVRYAjEmaZedXsPgWMBuKsxxmUwgDcDBmO12p2yyD4jINhF5VESaEjwXEblTRFpEpKW3t3fWQfePRpqwyu0udGOSdsmKBbhdwvN7Zv//pMmcbK87/hpYoqpnEallPJjoG6jq/ararKrN1dXVsw5oYDTahGUJxJhklRXksbahjI1tRzIdipmFTCaQDqApZrvRKTtGVY+o6oSz+T3g3HjPTZeBsQAikbtqjTHJu+C0Sra2DzAeCGU6FJOkTCaQTcAKEVkqIl7gVmBD7AEiUh+zeT3wuvP6KeBqEakQkQrgaqcs7QZH/ZQV5OF22TxYxszGBUsrCYSU196yfpBclbFRWKoaFJFPE/nhdwMPqOpOEfky0KKqG4C7ROR6IAj0AR91zu0Tkf9FJAkBfFlV++Yi7v7RgPV/GJMC5y6uRARe3neEi5ZVZTock4SMJRAAVX0ceHxS2ZdiXn8R+OI05z4APJDWAKcwMBagzEZgGTNrZQV5nFFXyiv75uRvP5MG2d6JnnUGR/1WAzEmRc5fWslrbw0QDIUzHYpJgiWQBA2MBaiwEVjGpMS6pjLGAiFae0cyHYpJgiWQBPUf9dtd6MakyFmN5QBsOziY4UhMMiyBJCAUVobGg7aUrTEpsrSqiBKfhy3tNhIrF1kCScDQmN1EaEwquVzCWY1lbLMEkpMsgSQgOpGizYNlTOqc1VjO7s5hu6EwB1kCSUB0Hqwyq4EYkzLrGssJhpVdnUOZDsUkyBJIAqIz8dowXmNSZ83CUgBetwSScyyBJGBgzJmJ15qwjEmZxooCSvI97O4cznQoJkGWQBIwYDUQY1JORFhVX8LuLquB5BpLIAkYdDrRS3wZnQHGmHlnVV0puzuHUdVMh2ISYAkkAcPjQQq9bluC05gUO6O+lOGJIO39Y5kOxSTAfgkTMDwesNqHMWmwqr4EsI70XGMJJAHD40FbSMqYNFhVV4II7O6yjvRcYgkkASMTQauBGJMGhV4PiyoL2dNtCSSXWAJJwJDVQIxJm+XVxbzZY7Py5hJLIAmwPhBj0mdZTTH7Dh8lFLaRWLnCEkgChseDlFoCMSYtllUXMREMc2jARmLliowmEBG5RkT2iEiriNw9xf7PicguEdkmIs+IyOKYfSER2eI8NsxFvJEaiDVhGZMOy6qLAWi1ZqyckbEEIiJu4D7gWmA1cJuIrJ502GtAs6qeBTwK3Buzb0xV1zuP69MdbyAUZjwQpiTfaiDGpEM0gbxpqxPmjEzWQM4HWlW1TVX9wMPADbEHqOpzqjrqbG4EGuc4xmOGx4OA3YVuTLpUFHmpLPJaAskhmUwgDcDBmO12p2w6dwBPxGz7RKRFRDaKyI3TnSQidzrHtfT29iYd7PB4dBoTa8IyJl0iI7GOZjoME6ec6EQXkQ8DzcDXYooXq2oz8OfAN0Rk2VTnqur9qtqsqs3V1dVJx2A1EGPSb1lNEa1WA8kZmUwgHUBTzHajU3YcEbkS+FvgelWdiJaraofz3AY8D5ydzmCHrAZiTNotqy6m76ifvqP+TIdi4pDJBLIJWCEiS0XEC9wKHDeaSkTOBr5LJHn0xJRXiEi+83oBcDGwK53BWg3EmPSLdqS3WS0kJ2QsgahqEPg08BTwOvBTVd0pIl8Wkeioqq8BxcDPJg3XPQNoEZGtwHPAPapqCcSYHHdadREAbb3WD5ILMvprqKqPA49PKvtSzOsrpznvJWBteqM7nnWiG5N+DeUFeFzCgT5LILkgJzrRs4HVQIxJP4/bRUNFAQeOjM58sMk4SyBxGh4P4MtzkWeLSRmTVosqC3mrzxJILrBfwzjZWiDGzI3FVYVWA8kRlkDiFEkg1nxlTLotrixicCzA4Ggg06GYGVgCidOQTaRozJxYVFUIYB3pOcASSJxsKndj5sbiaAKxZqysZwkkTraYlDFzY1FlJIFYR3r2swQSp+HxICX51oRlTLoVej1Ul+Rz4Ig1YWU7SyBxsk50Y+bO4kobiZULLIHEIRAKMxYIWSe6MXNkUZXdC5ILLIHEYcTuQjdmTi2uLKJraJzxQCjToZiTsAQSB5vGxJi5taiqAFU4NDCW6VDMSVgCiYOtBWLM3Gooj4zEau+3BJLNLIHEIVoDsftAjJkbDRUFAHRYDSSrWQKJg03lbszcqi3Jx+0SOqwGktUsgcQhWgMpthqIMXPC43ZRV+qzGkiWswQSh7drIJZAjJkrjRUFtPfbUN5sZgkkDjYKy5i511BRYE1YWS6jCURErhGRPSLSKiJ3T7E/X0Qecfa/LCJLYvZ90SnfIyLvSWecwxNBvB4X+R53Oi9jjInRWF5A19A4gVA406GYaWQsgYiIG7gPuBZYDdwmIqsnHXYH0K+qy4GvA191zl0N3AqsAa4Bvu28X1oMjwdsBJYxc6yhooCwQtfgeKZDMdOIK4GIyE0isldEBkVkSESGRWRoltc+H2hV1TZV9QMPAzdMOuYG4EHn9aPAFSIiTvnDqjqhqvuAVuf90sJWIzRm7tm9INkv3hrIvcD1qlqmqqWqWqKqpbO8dgNwMGa73Smb8hhVDQKDQFWc5wIgIneKSIuItPT29iYVqE2kaMzca7R7QbJevL+K3ar6elojSRNVvR+4H6C5uVmTeY+vvP9MxgPWDmvMXKov9wFYR3oWizeBtIjII8CvgIlooar+YhbX7gCaYrYbnbKpjmkXEQ9QBhyJ89yUaawoTNdbG2Omke9xU1OST8eADeXNVvE2YZUCo8DVwPucx3WzvPYmYIWILBURL5FO8Q2TjtkA3O68vhl4VlXVKb/VGaW1FFgBvDLLeIwxWaahosD6QLJYXDUQVf1Yqi+sqkER+TTwFOAGHlDVnSLyZaBFVTcA3wf+TURagT4iSQbnuJ8Cu4Ag8ClVtXmfjZlnGsoL2NExmOkwzDTiSiAi0gh8C7jYKXoR+EtVbZ/NxVX1ceDxSWVfink9DnxwmnO/AnxlNtc3xmS3heUFPL2rG1UlMgDTZJN4m7B+QKTZaKHz+LVTZowxaVNX6mMiGKZ/NJDpUMwU4k0g1ar6A1UNOo8fAtVpjMsYY1jojMTqHLR+kGwUbwI5IiIfFhG38/gwkdFQxhiTNnVlkXtB7G707BRvAvk4cAvQBXQSGRGV8o51Y4yJVV8WqYEcsgSSleIdhXUAuD7NsRhjzHEWFOfjcQld1oSVlU6aQETkb1T1XhH5FnDCXdyqelfaIjPGnPLcLqG21Een1UCy0kw1kOj0JS3pDsQYY6ZSX+ajc8ASSDY6aQJR1V87L0dV9Wex+0RkyvszjDEmlerKfOw8NNvJv006xNuJ/sU4y4wxJqXqy3wcGhgjMouRySYz9YFcC/wZ0CAi34zZVUpkChFjjEmr+rICJoJhBkYDVBR5Mx2OiTFTH8ghIv0f1wObY8qHgc+mKyhjjImKDuXtHBy3BJJlZuoD2QpsFZGfOAs6GWPMnKore/tu9NULZ7uOnUmlmZqwfqqqtwCviUhsA6QAqqpnpTU6Y8wpb2F55G50G8qbfWZqwvpL53m2a38YY0xS3r6Z0BJItjnpKCxV7XReHgYOOnek5wPriPSPGGNMWkVvJjxkd6NnnXiH8b4A+ESkAfgt8BfAD9MVlDHGxKor81kNJAvFm0BEVUeBm4Bvq+oHgTXpC8sYY95mCSQ7xZ1AROQi4EPAb5wyd3pCMsaY4y0sizRh2c2E2SXeBPIZInee/9JZj/w04Ln0hWWMMW+rKytgPBC5mdBkj7gSiKr+XlWvB+4TkWJVbZvNTLwiUikiT4vIXue5Yopj1ovIn0Rkp4hsE5H/FLPvhyKyT0S2OI/1ycZijMl+0ZsJu4asGSubxJVARGStiLwG7AR2ichmEZlNH8jdwDOqugJ4xtmebBT4iKquAa4BviEi5TH7/1pV1zuPLbOIxRiT5WpLLYFko3ibsL4LfE5VF6vqIuCvgH+dxXVvAB50Xj8I3Dj5AFV9Q1X3Oq8PAT3YOuzGnJKO1UCsIz2rxJtAilT1WJ+Hqj4PFM3iurUx95h0AbUnO1hEzge8wJsxxV9xmra+LiL5Jzn3ThFpEZGW3t7eWYRsjMmU6pJ8RCyBZJt4E0ibiPydiCxxHv8DaDvZCSLyOxHZMcXjhtjjNDKsYtqhFSJSD/wb8DFVDTvFXwRWAecBlcAXpjtfVe9X1WZVba6utgqMMbkoz+2iujjfEkiWiWtNdODjwD8AvyDyY/+iUzYtVb1yun0i0i0i9ara6SSInmmOKyUybPhvVXVjzHtHay8TIvID4PNxfg5jTI6qK/NZH0iWmWkyRR/wSWA5sB34K1VNxTi6DcDtwD3O82NTXNsL/BL4kao+OmlfNPkIkf6THSmIyRiTxepKfRw4MprpMEyMmZqwHgSaiSSPa4Gvpei69wBXiche4EpnGxFpFpHvOcfcArwT+OgUw3V/IiLbnbgWAP87RXEZY7KU1UCyz0xNWKtVdS2AiHwfeCUVF1XVI8AVU5S3AJ9wXv8Y+PE051+eijiMMbmjttTH4FiAMX+IAq9NhJENZqqBHGuusgWljDGZZDcTZp+ZaiDrRGTIeS1AgbMdXVDKlgczxsyJutK3VyZcumA2dxGYVJlpSVurJxpjskJ0adtuq4FkjXjvAzHGmIyqO3Y3+kSGIzFRlkCMMTmh0Ouh1Oehy1YmzBqWQIwxOcOG8mYXSyDGmJxRV1Zg05lkEUsgxpicUVeabzWQLGIJxBiTM+rKCugdniAYCs98sEk7SyDGmJxRV+ojrNA7YiOxsoElEGNMzojejd5p/SBZwRKIMSZnRJe27bYEkhUsgRhjcobVQLKLJRBjTM4oL8zD63HZdCZZwhKIMSZniAj1ZT6rgWQJSyDGmJxSW2p3o2cLSyDGmJxSX+azu9GzhCUQY0xOqXNqIKqa6VBOeRlJICJSKSJPi8he57limuNCMeuhb4gpXyoiL4tIq4g8IiLeuYveGJNJtaU+/MEw/aOBmQ82aZWpGsjdwDOqugJ4xtmeypiqrnce18eUfxX4uqouB/qBO9IbrjEmWxxb2jaNzViHBsa477lWXm47krZrzAeZSiA3AA86rx8Eboz3RBER4HLg0WTON8bkttpja6OnZ12Q/qN+bvr2S3ztqT3c+q8beeb17rRcZz7IVAKpVdVO53UXUDvNcT4RaRGRjSISTRJVwICqBp3tdqBhuguJyJ3Oe7T09vamJHhjTObUp3llwv/3wpt0D4/zyJ0XcnptCX/3qx34gzZ541TSlkBE5HcismOKxw2xx2mkJ2y63rDFqtoM/DnwDRFZlmgcqnq/qjaranN1dXXiH8QYk1Wqi/NxCWlZmXBkIsiPXjrA9esWcsFpVfz1e07n0OA4z+62WshUPOl6Y1W9crp9ItItIvWq2iki9UDPNO/R4Ty3icjzwNnAz4FyEfE4tZBGoCPlH8AYk5U8bhfVJelZF+SZ17sZC4T4iwsXA3DZ6TXUl/l46JWDXHNmfcqvl+sy1YS1AbjdeX078NjkA0SkQkTyndcLgIuBXU6N5Tng5pOdb4yZv+pK03M3+m+2dVJX6uOcRZGBoW6XcNM5Dby4t5dBG/V1gkwlkHuAq0RkL3Cls42INIvI95xjzgBaRGQrkYRxj6rucvZ9AficiLQS6RP5/pxGb4zJqLoyX8rnw/IHw7y49zBXr6nF5ZJj5VecUUtY4fd7rQ91srQ1YZ2Mqh4BrpiivAX4hPP6JWDtNOe3AeenM0ZjTPaqK/Xx0pupHWK7tX2AsUCIdyxbcFz5usZyKou8PL+7h+vXLUzpNXOd3YlujMk5dWUFDI8HGfUHZz44Tn968wgicOFplceVu13CJcsX8IfWw3b3+ySWQIwxOaeuLB9I7c2Er+zr4/TaEsoLT5zY4ryllfQMT9Den5qRX+Gwsr19kN7h3F6a1xKIMSbn1JUWAKlLIKrK1vYBzlk85axKnLckUr5pf9+srxUMhfnPP97M+/7lD1x677O8mMN9K5ZAjDE5p+7Y3eipSSAHjowyPB5kbUPZlPtX1pRQ4vOwaX//rK/10KaDPL2rm/92+XKWVBXx2Ue2pLQpbi5ZAjHG5Jy60tQubbu9YxBg2gTicgnnLq5g84HZ1UDGAyHue7aV5sUVfO6qlXzl/Ws5POLnwZcOzOp9M8USiDEm5xR43ZQV5KVsKO+OjkG8Hhcra0umPaZ5cQVvdI8wMOpP+jpP7+qma2icu65YgUgkKV26YgEPvrSfUHh2HfT+YJjNB/rpP5p8fImyBGKMyUmpXNp2W/sgZ9SV4PVM/5MYvblwa/tg0td5bEsHdaU+Lln+9lDhW89bRNfQOH+axbDk4fEAN973Rz7wnZd4573Pzeq9EmEJxBiTk2pLU3MzYTis7Dg0yNrGqZuvos5sLEMEth4cSOo6Q+MBfv9GL+9bVz/pRsUaSnwefrUl+RmZ/vsvd7Cne5gv37CGmtJ87nr4NYbH03/nvCUQY0xOSlUNpGNgjOHxIGsWnjyBlPryWFZdnHQC2bSvj0BIueKM4ycf9+W5effpNTy3u4dwEs1YOzoG+fXWQ3z63cv5yEVL+Kdb1tM7PMF9z72ZVJyJsARijMlJtaU+Do9MEAjNbqr1vT3DAKysLZ7x2HWN5WxtH0jqhsKNbUfwelysbyo/Yd8VZ9Rw5Kifre2JJ6fvPP8mJT4Pd1y6FID1TeW8b91CfrLxACMT6R3dZQnEGJOT6sp8qELPLG/Ge6N7BIDlNdN3oEetbyrj8IifjoHEbyh8eV8fZzeV48tzn7DvXSurcQk8t3vKicmn1XfUz293dXFLcxOlvrxj5R+/eAnDE0F+8Wp7wnEmwhKIMSYn1aVoads3uoepLc2nrCBvxmPXObWHrQcT60gfGg+wo2OQC06rmnJ/eaGXcxdX8OyexBLIhi0dBPcexI8AABE2SURBVELKzec2Hld+9qIKVteX8vPNlkCMMeYE0XtBZptA9naPnHT4bqxVdaV43a6Em5pa9vcR1hPn2Yp12ek17OgYSmh6k99s72RVXQln1JeesO/GsxeytX2QfYePJhRrIiyBGGNyUn0K7kYPh5XWnhFWxNF8BeD1uFi9sJQtCXakv9zWh9ftOjYUeCrvXBFZMfWPrYfjes8jIxNsPtDP1Wvqptx//boGRCJDh9PFEogxJieVFeSR73HNamnbjoExxgIhVsTRgR61vqmc7e2DBBPovN/YdoT10/R/RK1ZWEpFYR4vxDk31jO7ewgrXL26dsr9dWU+LlxaxYYth9I2i7AlEGNMThIR6st8dA0l34n+Rnf8I7Ci1jeVMxYI0do7Etfxw+MBdhwa4oKTNF9BZLqUS1ZU8+Le+KaNf3pXNwvLfKxZeGLzVdT71i2k7fBRXu8cjivWRFkCMcbkrLoyH51JjIiKSmQEVlS0I33LW/E1Y7Uc6CcUVi6cpgM91qUrFtA7PMGe7pP/4I/5Q7y4t5crV9ciItMe9541tbhdwm+2H4or1kRZAjHG5KyG8sKkhtRG7e2JfwRW1JKqQkp9nrg70je2HSHPLSft/4iK9oO88MbJm7H+0HqY8UCYq6ZpvoqqKs7notOqeHx7V1qasTKSQESkUkSeFpG9zvMJ/7Ii8m4R2RLzGBeRG519PxSRfTH71s/9pzDGZFpTZQFdQ+NMBENJnZ/ICKwoEWFdUzlb4hzK+3JbH+sayynwTt//EVVX5mNlbTEv7j15R/rTu7ooyfdwwdKZazXvPauefYePsvPQUFzxJiJTNZC7gWdUdQXwjLN9HFV9TlXXq+p64HJgFPhtzCF/Hd2vqlvmJGpjTFZprChEFQ4NJD4SK9ERWLHWN5XzRvfwjOt4jEwE2d4xGFfzVdSlK6p5eV8f44Gpk2IorDzzeg+Xrao56eSPUdeeWccPP3Yep9cl/jlnkqkEcgPwoPP6QeDGGY6/GXhCVUfTGpUxJqc0VURWJjzYl/hPQ3t/ZARWIh3oUesaywmFdca/6lv29xEK64wd6LEuXbEAfzDMK/umXnvklX19HDnq59ozpx6+O1l5oZfLTq8hz536n/tMJZBaVe10XncBJ2/Ig1uBhyaVfUVEtonI10Ukf7oTReROEWkRkZbe3txdOtIYc6KmykKApNYqj86BlcgQ3qizmiITL840seLL+/rwOItRxeuCpVV43a5p+0Ge2NGJL8/FZadXxx9wmqQtgYjI70RkxxSPG2KP00jPzrS9OyJSD6wFnoop/iKwCjgPqAS+MN35qnq/qjaranN1deb/wY0xqVNb6iPPLRzsT7wGsrcn8RFYUTUlPhrKC2a8oXBj2xHOaiyj0OuJ+70LvG4uWlbFkztP7PgOh5Und3Rx2cqahN4zXdKWQFT1SlU9c4rHY0C3kxiiCeJkE8DcAvxSVY9Nbq+qnRoxAfwAOD9dn8MYk73cLmFheUFSTVh7u0eoKUlsBFasdU1lJx2JNTIRZHv79PNfncx1Z9XT3j92wuJVr+zvo2d4gmvXxtd8lW6ZasLaANzuvL4deOwkx97GpOarmOQjRPpPdqQhRmNMDmiqKORgEk1YrT3DSTVfRa1rLOdg3xhHRqa+kXHjm0cIhpVLY1YfjNfVa+rwul38ctJsug+98hYlPg9Xrz61E8g9wFUishe40tlGRJpF5HvRg0RkCdAE/H7S+T8Rke3AdmAB8L/nIGZjTBZqqiygPcEaiGryI7CiojcUbptmids/tB7Gl+fi3CXx939ElRXkcd1Z9fxsczuDo5HGl+6hcZ7Y3sUHzmmMa0jwXMhIAlHVI6p6haqucJq6+pzyFlX9RMxx+1W1QVXDk86/XFXXOk1iH1bV+OYUMMbMO40VhRw56p9xSG2szsFxjvpDLK9JvgaytqEMl8Crb/VPuf/Fvb1csLSKfE9yP/Z3XLqUUX+Ibz/fCsC9T+4B4GMXL0nq/dLB7kQ3xuS0RmcobyIjsd7uQE8+gRTlezirsZwXprjp79DAGG/2HuXSFYk3X0WtWVjGbecv4v4X2/gvP97Mz19t545Ll7K4qijp90w1SyDGmJwWHcqbSEf6XmeuqRWzSCAAl6+qYVv7AIcn9YM8vycyBPeSWSQQgC9dt5qrV9fy/J5ePnhuI5+7auWs3i/VMj8OzBhjZqGpIvEE0tozQmWRl6riaW8hi8vlq2r456ff4NndPdzS3HSsfMPWDk5bUMTpCU6TMlmB1813/6J5Vu+RTlYDMcbktAXFXgry3AmNxNrbMzKr5quoNQtLWVJVeNzSsV2D47y8r4/3rVt40ply5wNLIMaYnCYiLK4qjHvp1rdHYM0+gYgIt5zXxMv7+mh1+lUe+OM+BLjpnIZZv3+2swRijMl5y2qKeTPOBZ56RyYYHAukpAYC8MFzmyjyuvmHX+/kYN8oP/rTfm5c35BVnd3pYgnEGJPzllUXc7BvNK5p3VudRaRmcw9IrOqSfO6+dhUv7j3Mu772HHluF5/Nss7udLFOdGNMzltWXURY4cCR0RnX94gO4Z3NXeiTffjCxeR73GxtH+AvLlp8bGTYfGcJxBiT85ZVR5LBmz0zLxC1u2uYUp+HmpLZjcCKFe0LueW8ppkPnkesCcsYk/OWLoj0N8TTD7Krc4jVC0vn/QipuWAJxBiT84ryPSws8/Fm78lHYgVDYXZ3DrFmYdkcRTa/WQIxxswLp1UX0zZDDWT/kaNMBMOsri+do6jmN0sgxph5YXlNMXt7RgiHp12f7tgStKsXWgJJBUsgxph5YXV9KaP+EAdOMqXJrs4hvG7XsU53MzuWQIwx80K0VrHz0NTrcwDsOjTEitpivB776UsF+1c0xswLK2qL8biEXU4z1WThsLL14ABnNVoHeqpYAjHGzAv5HjfLa4qP9XNM1nZ4hKHxIGcvSnyFQDM1SyDGmHljfVM5Ww4OTNmRvvlAZOXAcxdbAkmVjCQQEfmgiOwUkbCITDvZvYhcIyJ7RKRVRO6OKV8qIi875Y+IiHduIjfGZLPmJZUMjgVonWI47+YD/ZQX5nHagvk/yeFcyVQNZAdwE/DCdAeIiBu4D7gWWA3cJiKrnd1fBb6uqsuBfuCO9IZrjMkFzU7tYtP+vuPKVZU/th7hgqWVdgd6CmUkgajq66q6Z4bDzgdaVbVNVf3Aw8ANEvn2LwcedY57ELgxfdEaY3LF4qpCqkvy2dh2fAJ5s3eEjoEx3rWyJkORzU/Z3AfSAByM2W53yqqAAVUNTiqfkojcKSItItLS29ubtmCNMZknIly2sprn9/TgD4aPlUfXKH/nytmtUW6Ol7YEIiK/E5EdUzxuSNc1p6Kq96tqs6o2V1dXz+WljTEZ8J41dQyPB9nYduRY2WNbDnFGfSmNFafGNOtzJW3TuavqlbN8iw4gdm7kRqfsCFAuIh6nFhItN8YYLlmxgJJ8Dw9veot3rqxm16EhtncM8vfvWz3zySYh2dyEtQlY4Yy48gK3AhtUVYHngJud424HHstQjMaYLOPLc/ORdyzmiR1d7O4a4pvP7MWX5+LGs+f/GuVzLVPDeN8vIu3ARcBvROQpp3yhiDwO4NQuPg08BbwO/FRVdzpv8QXgcyLSSqRP5Ptz/RmMMdnrjktOo6LQy3Xf/ANP7uziritWUF5oo/1TTSJ/0J8ampubtaWlJdNhGGPmwO6uIe5/oY3V9aV8/OKluFw2fDdZIrJZVU+4Z8+WtDXGzEur6kr551vWZzqMeS2b+0CMMcZkMUsgxhhjkmIJxBhjTFIsgRhjjEmKJRBjjDFJsQRijDEmKZZAjDHGJMUSiDHGmKScUneii0gvcCDJ0xcAh1MYTi6wz3xqsM88/8328y5W1ROmMz+lEshsiEjLVLfyz2f2mU8N9pnnv3R9XmvCMsYYkxRLIMYYY5JiCSR+92c6gAywz3xqsM88/6Xl81ofiDHGmKRYDcQYY0xSLIEYY4xJiiWQOIjINSKyR0RaReTuTMeTaiLSJCLPicguEdkpIn/plFeKyNMistd5rsh0rKkmIm4ReU1E/sPZXioiLzvf9SMiMq/WQRWRchF5VER2i8jrInLRfP+eReSzzn/XO0TkIRHxzbfvWUQeEJEeEdkRUzbl9yoR33Q++zYROSfZ61oCmYGIuIH7gGuB1cBtIrI6s1GlXBD4K1VdDVwIfMr5jHcDz6jqCuAZZ3u++Uvg9ZjtrwJfV9XlQD9wR0aiSp//CzypqquAdUQ++7z9nkWkAbgLaFbVMwE3cCvz73v+IXDNpLLpvtdrgRXO407gO8le1BLIzM4HWlW1TVX9wMPADRmOKaVUtVNVX3VeDxP5UWkg8jkfdA57ELgxMxGmh4g0Au8FvudsC3A58KhzyLz6zCJSBrwT+D6AqvpVdYB5/j0TWbq7QEQ8QCHQyTz7nlX1BaBvUvF03+sNwI80YiNQLiL1yVzXEsjMGoCDMdvtTtm8JCJLgLOBl4FaVe10dnUBtRkKK12+AfwNEHa2q4ABVQ062/Ptu14K9AI/cJrtviciRczj71lVO4B/BN4ikjgGgc3M7+85arrvNWW/aZZAzDEiUgz8HPiMqg7F7tPIeO95M+ZbRK4DelR1c6ZjmUMe4BzgO6p6NnCUSc1V8/B7riDyF/dSYCFQxIlNPfNeur5XSyAz6wCaYrYbnbJ5RUTyiCSPn6jqL5zi7mjV1nnuyVR8aXAxcL2I7CfSLHk5kf6BcqepA+bfd90OtKvqy872o0QSynz+nq8E9qlqr6oGgF8Q+e7n8/ccNd33mrLfNEsgM9sErHBGbXiJdMBtyHBMKeW0/X8feF1V/zlm1wbgduf17cBjcx1buqjqF1W1UVWXEPlOn1XVDwHPATc7h823z9wFHBSR052iK4BdzOPvmUjT1YUiUuj8dx79zPP2e44x3fe6AfiIMxrrQmAwpqkrIXYnehxE5M+ItJe7gQdU9SsZDimlROQS4EVgO2/3B/x3Iv0gPwUWEZkG/xZVndxRl/NE5DLg86p6nYicRqRGUgm8BnxYVScyGV8qich6IoMGvEAb8DEif0jO2+9ZRP4B+E9ERhu+BnyCSJv/vPmeReQh4DIi07Z3A/8T+BVTfK9OIv0XIk15o8DHVLUlqetaAjHGGJMMa8IyxhiTFEsgxhhjkmIJxBhjTFIsgRhjjEmKJRBjjDFJsQRiTBJEpEpEtjiPLhHpcF6PiMi303TNz4jIR06y/zoR+XI6rm3MVGwYrzGzJCJ/D4yo6j+m8Roe4FXgnJg5nCYfI84xF6vqaLpiMSbKaiDGpJCIXBaztsjfi8iDIvKiiBwQkZtE5F4R2S4iTzrTxyAi54rI70Vks4g8Nc3MqJcDr0aTh4jcJZH1W7aJyMNwbL6j54Hr5uTDmlOeJRBj0msZkR//64EfA8+p6lpgDHivk0S+BdysqucCDwBTzXRwMZFZZKPuBs5W1bOAT8aUtwCXpvxTGDMFz8yHGGNm4QlVDYjIdiJT4TzplG8HlgCnA2cCT0daoHATmXZ8snqOX/hqG/ATEfkVkSkronqIzDprTNpZAjEmvSYAVDUsIgF9u9MxTOT/PwF2qupFM7zPGOCL2X4vkcWh3gf8rYisdZq3fM6xxqSdNWEZk1l7gGoRuQgi0+qLyJopjnsdWO4c4wKaVPU54AtAGVDsHLcS2DHF+caknCUQYzLIWSb5ZuCrIrIV2AK8Y4pDnyBS44BIM9ePnWax14BvOkvTArwb+E16ozYmwobxGpMjROSXwN+o6t5p9tcC/66qV8xtZOZUZQnEmBzhLARVq6ovTLP/PCCgqlvmNjJzqrIEYowxJinWB2KMMSYplkCMMcYkxRKIMcaYpFgCMcYYkxRLIMYYY5Ly/wGvVuX8GIreuQAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2de5yV5XXvv4tBRgYZYRiEASGMjoqKCgHFaCR4jZeYxB6T2NxzYk2a5vTeU9uc06Y9pz323rTxJLG5NidWjYkpNqYaL0SrCREq8YoKQkRBEFAHGDLI8Jw/9kYHnIHnC3vP8ML7+3z8OOxZ8+x3v3vvZz1rrd/6rUgpUaJEiRIlSvSHIYN9ASVKlChRYv9G6ShKlChRosRuUTqKEiVKlCixW5SOokSJEiVK7BaloyhRokSJErtF6ShKlChRosRuUTdHERFfjYi1EfFor8daIuKHEfF09f+jq49HRPxDRCyNiIcj4s31uq4SJUqUKOFQz4ji68CFuzx2NXBXSukY4K7qvwEuAo6p/ncV8IU6XleJEiVKlBCom6NIKd0LbNjl4XcB36j+/A3g3b0e/+dUwU+AURHRVq9rK1GiRIkS+Rg6wM83LqW0uvrzC8C46s8TgZW97J6rPraaXRARV1GJOmgawcxjpub5um6GqQvtYoSyH8nGbNtDX9mq1t52eL4/f5pj1NpjeVHZt67d1ffvBg1qaV4ac7j7A4FXcGu/yiHK/pjuZfnGnWppeFXa/0LY2qPieGG7Ta7tbjkvNzYr+2WLOtellMa6Z9kZHRGpK9N2NdyeUto1q1JIDLSjeA0ppRQRWj8kpXQdcB3AcbMOS3+38OSsvxvOFvU8PXKXm7STn9s9Ji9dq9Ze2jEh2/bX+XO19vu4Udl/5C9vyjc+VC3NX/36h5X9MLqzbZfRodbeyEhl/7Wlv5pv/EW1NDRKe+OIVrmlt30t33boE25t6cv5xlS3B380bvq5e4Y3ogv4RKbtZ6F1X59vf8FAO4o1EdGWUlpdTS3t2DGfByb1sjuy+thusZGR3MPZWU98NveoCz17+Y+V/br2w7JtH+44Vq39MqOybWewWK09mpeVvfrEyBPlCqYo+23CmU94Y3C6W5zA48peBaDuLYJp0t4ElU+7pV9uzv+ctx6+yS0u0TMI59xgEE/Xg4iBfs3zgI8A11T//6+9Hv90RNwAzAZe6ZWi6hebGcFCZmY98Rk84K70bmc+8pfzvxQPNJ2h1h5ObrALZ3GfWnsM65Q9Rzhzgy6alL1xLNaB2vuYxKk/nlVLQ7u0zw9AwWVvFJL8rMRmZ99qP7s1wBBg+IA/6+Cjbo4iIv4FmAu0RsRzwB9TcRA3RcTHgZ8D762a3wZcDCylEt19LOc5tnIIK3cKRHZn62oUPOLMG9+Zb/tQ03S1tjkNX9ElUkNAd6O8L+a0KlNPTcIhAtnvPcAoGTmduETUHMClTXrc0iLDVoGJbo53S4/qzD8QvdzittSGw92NmcIKZV8LBLqUckCgbo4ipfTL/fzq3D5sE/Br9jkqYWDeh6vBfjvl6cYUEFfII2LuawRo/IlamsYxsrAuNpaVzY64Nob1yr5bJO+to+B+Z67SQzYqM8VpcCm/093SQ0U0tH7aGLe4JD8c3b3U/UENUKaeCohDeJW2zGqcKXwCID/j5k52yeDVbIg252yLmUtPeVO27QJmq7XtZj6Kl+q2Ns84cwWbSrLsIXHISRe4pWNhvm3XNJdKtKnHjuXPKftaoJYRRURcCHyOiov8ckrpml1+/0kqB+geYBNwVUpJFs9qg0I7isN5hbdze5at3ijkF2hzSz7PcByO9WQ2RJ2mWO7Mn+K4bNvbuFitfR53KfvZ/DTbVt1DeGMH0J5gHPRb5doPSftX8k0XtrhK+ambH92zURWWOWhIG4D/rNcAtYooIqIBuBY4n0o7wIMRMW8XR3B9SumLVft3An/LG5uYBwSFdhTNdGY7ikZciuWFcxxXz3zIj+NJtbah3uq0hiysrnmt9WXPeJwT1NpnM1/Zn8W92batMq2lvxmmVv4hubaM+kwEsknSgE2kPUx+5zTqSKzoDzWMKE4DlqaUngGoEnneBa/T7VJKvYnOI4BBG0daaEdxaNdWpi7Ko0a/NNOle+xp2OTXZyLid2QEYtMaMqIwabP1kkZuexfezS3Ztq2rJVXT9VvSKXx5s1xb1OwrEOvb/iJDZrDkBNO0CvBsm/UULpLvCzVkPfXVZPyGXG1E/Brw28Aw4JzaPLVHoR0FnZCbrRjV7r4QD7Q4CuuxIkq4nO+otcd15juh5BqzCZGmgPpy162jGH+3uHh72XIzP8RkWWzK5ChpL6CiVWBJW36NykZxNj1sI9ZaOAoZUbRG7FTVua7aMJyNlNK1wLUR8X7gf1BpKxhwFNtR/ALI7P60HO0nW/Jz8eDC7ImPyjSI+C4/e5E7ZY28wJ3iDCbInImmMJvAzKqDyG7o4aZ3waaSzpT2IvU0cbX7LM5ry+tbApjUfadae1SjcxT3cZayrxXEprkupTSrn9/ZJuMbGESx1GI7iu14GmsmbGFtvUneWmaSsF9/kaNrPYCLnAyOxtEXdU7bbLju4OzTPaY1RtaFuvsjmveDBuEorMzGmrb8GpWFTVXZCLQWqGGN4kHgmIhop+IgrgDev9NzRRyTUtrx7b8Ev3PUDMV2FEMhu7Yqm7/spqVEB20ELDY5SzG8nbcr+5NEJ6Kxhb2QEzGtMVZ3yKa/TT+C/LpvbMqXzQD3GZj8ivswmqhva6MLyxp6XK9TU4NzLLVArVhPKaVtEfFp4HYq9NivppQei4g/BRamlHaoVZxHRRbyJQYp7QRFdxSHk01jXTfWfdms1o8qxFmRN3FCtOkb6yhMlGBZTJrCbNI989zSWl9JdOZbR2E0rQDWCC83eYxzFKamYa4DYGSDIxyYumCtUEsJj5TSbVRUKXo/9ke9fv6NGj3VPqPQjmLb4UNYd2ne6WmV2lVglmQmqc5vdyneXmD1Krf4pgn54f5JPS6iaH5app5MesimKKV9EtcS8v0c2e3qSD2N+V/rbinhYTbntYJKDZ4ocYIOE/cdpYRHAfEKh/N9zs+ytU1us1mg7JWjkMwks7Tq4gaGHeooOCat0fyE3PjvcObqPtpitmQmbRWpzcb+ypv9YMSG7cq+sS3/4h9vcp5i+ob8zfkpSQixMjsnrpR6XDVAKeFRQGxkJPMzZcateuyFq+e7ixF6PEaSHGDkuPyQ3MqDTGpxVV5Vi7EFZKuvZPY4W5yWukONgqm7TbKYjL4SwPC2/Nz9I+TNc9mBGSvzHUV3izu02M+u7p6vAcqIooDooonFmXQTy8DRMwPW5JuubHe71pqm/BDeMkFs85/qcJc9GjqTYJidJ8m1LQQD65GZ7hQ/o9vdmJHkHyzW2qq9OBA1yrBM9+jYWl8NUEYUBcSrHJJde9BSBVYUThQojQwGwHzmZttaBVYbaW0xrCqremrZYKYo7JTdNYXVOIqlM920vRmSw2poprZQbljgttPaXsu2Otbu+kMZURQQiaBne96Hq3uIPH7YsFZRWF2Ifb/odbC1lQ9wvbI3MyA05J6lIhDDStobCP+sUyzy0G9opq0N7mCxriM/bXqE9Py2d+nJ5qOVPex7TSMoBxcVDsPZwvFD8mismnppi59Wv0dgUWd+9bO9eYVa2+ScQZ5uZe+KlqoQ6cEkNbDCmavoxhIrlrYdqeyP6Mlf3/a6mAbNuT3z1dr2oJA73fJ11MZRHJK7a1p5+P0YhXYUI9icfYK2mjb64CzolGYQEcCWTfnpnq5m13CnUyz5KtPeeVo2mL12Ae1YRD2mQ9bLluJSVYc15Kd8pvSsUGvf2XBetu2lrzga28gWl6qy9NtaIAKGlo6iWBjJxuymLpsvtSMi6ym/fNio/GtvsJ9O219gvvtOgNcr34qoL2Qq8YUOF1KOPybfU3S86Abu3D/WyayYOtL4Zx3j4OV2OTNCQE+hHARESAHIAwSFdhTNPRt5e2eefOzLzY6S+sJUt1H0TK3frexoyj+B6hSbhRm16gbceUdhIhAZfazvcBLpPZPy3/+J97u6wKaxjoihmG/d+fPYwRXK9cxsecixxI1aQEUUBxAK/ZLjFRiaKc3QepKTB7jnlLco+2dFrspGN6eJArVOsck6wgrB7pliVVIthdXoK9Vx7ga4QuxEIRsPfmaEchQyPTiFFdm2Vg3Bfi/sALBaIAIOGQRa7mCj0I6CTeQ3acmMzJOnuK5SVeSTGkhGM0nPBq/nh15uzrqmMUfY5g/DA3wdabXJPcrXaeXajZNbN8lF2qb4bRlyU+QHZjAcxcHaSFHsl7yFfIpki1t6BVOU/Z3kF/mOluyLz3T+RbZtd2P+7G5AO9Apk4WxdRTyWp6fk785T1xQ3yE6T4pZ4t1SwmO6HJq9iPwnsMN/5izNn1O+oMPlHsewTtlPXZI33bKmKB1FAZHI31zkQdtqJq1+Nv/0tH6ymxkxVNQFhvY4XSANs8m5pm/9HpkT68TJzlEc0S1a7YGmxvzc/cNNTpr21JWGagYNk/KjId0XsyrfUWztcErGqpkT/KTAWqHYu+ZeodgveRj5yqqyL0Kzhzbl30o9cMWczO1oEyunbRyFlUGRDCwlqS4p9yPWOoc7ZlK+I7Kn+FNXOUfROCl/B63n8B9bc+ixjRTNzrwmCHxj6AGAYjuKJmBGpq08OBm9HABG5X859RfCUDutsJ5tLDSOwvaiyFtuaMZL2vNnPQNM/ZlLa4yZlJ82uUdIsuwNDhMbtKakivfU1lbstVhxTf0B6wtl6qmAaCZ7cJHdEFtlvvToCXWksIosWKc8xTdPdfa8R9hafSWZSjAbkR35OnWbcxStgqqpxe+kppGZFGg/i8vb27JtDUMKPEvKRmaQnzbrF8GgiBEONgrtKLqbDmHpzLzuTFtzsBxto8KqKazCyT0vhfiarRS4Yey8Va4tafETOvOlKlY2y/BGHixGduWfVkc1vaTWfnaSE3sym7/9LD4kvP9lG36g1l7X4mp3ViK9Zo6iRrtmRFwIfI5KMuvLKaVrdvn9bwNXUqnEvgj815TSIFTwC+4o1jOGb/GBLFsreWzF9d7Hjdm25vQJqHBfv6G210E08r50uetFGH2D6xcYKmo33ae4wqrtzB63If/GzGhyYZ+Nhi7r+l62rWnmBDc695fWOkfR1OLefy2RXgvUyFFERANwLXA+8BzwYETMSyn1Fq97CJiVUuqKiF8F/hJ4374/u0ehHcVmRrCA07Js1RwF4CzuU/bn9eR1iAM0bJOFcnGKP8YdyjxzRDiWx6e51MCZExa5axFOywrxLZP6Sq1D86/dzmO/l7OUvcGoTpe339IsnL8eL+E+jMNFl3hNUZti9mnA0pTSMwARcQPwLnj9w5FSuqeX/U+AD9bkmfcChXYUWxjO4sxqtpEeAF9Ya54nHJFVVTU9IPWeuyAyFTbnzDHSUYiZIbY5aylOwvqE5vzNv3W525w3tTtmUldT/mY++kV3ijdprc2TXE+P/c5NwMmP1AQuomiNiN456etSStdVf57Izt+m59i96M3HARei1RCFdhSvvjqM1avyNqMhDe4U3zNOHht+JGztXTfienLMph4uJCIKy+5a3pZfKAVoX5+/UdRTThuchMfotW5zHt7u7NeLLvGRjU6g8FjhcJc1uqisCBIe0lGsSynJ9so+njLig1T4hm/b17X2FoV2FPQEbMorUm+Xxeyt41xOW8lv283ZRAnWUdi+C9HrMM7Mh8Wf4qdMyHcUo3/mNtuuU1zz1zqxObcf6k7CVpbcRHKHNbvNeYZojlkgVSFtXfC4rqeUfU1QO9bT8+xcfTyy+tjOTxdxHvAZ4G0ppcFqMSy2oxgybBtNR+bRWLs2usKqHhEpTtpr5GzocYI9tO0P3dpDraMQPR32FH+jrNO1t6zItu1Y4E7O405xTs7MRtjW7hrobL3sPlHTsHn+U3+Wf+03njJFrW2lShrtZ7cWqB3r6UHgmIhop+IgrgDev9NTRcwAvgRcmFKyg4JrikI7isYh3dmsjZea6qejX7mYfFN3PoRxIhe/rlnOUTjdzSMwshytj7pc/MZpLhdvTs4ddUyxAKwnX5Z8abNs/lvqGJEbO/Lvo7luABMk2sPW1iI0KNTIUaSUtkXEp4HbqZTHv5pSeiwi/hRYmFKaB/wVcBjw7YgAeDalVO+hvn2i0I5iBJuZSV4B1FLpmqS0sxm2JVf2FFaBl2ZKCut6cfXyxDdsmmOmmXSPravbaOg2Lsm2XSEHb0x9RlLnRWlAz+8WRAxLA+82kiyg39OaoIYSHiml24Dbdnnsj3r9nK80WmcMiqOIiN+i0kiSgEeAjwFtwA1UZsUtAj6UUtrtzjGCzZzBA1nPaQp8sBfd02KIzvDb3dIm3TO6y0UIC5ucCNKZxwtmkpuEqSUfXmZ0tq0dRDV+ibuPDVPzyRJaiE9mpk2zqO4SFz7OFpttRLF8rCM/UAuWVCnhMTCIiInArwMnpJS2RMRNVPJzFwN/l1K6ISK+SIUO9oXdrVUZhXrP7kxegxU/s5sWQgrDcUFQ8tuNIk0FsGLaFGU/ZdKKbNuJje5EaTcWUyx9Qso9jF/+Y2XfNDU/0jIMKUCfnE23tY20l046MtvWRmWPyMlVRtq9gho5igJkyGqNwfKNQ4HhEfEqFWm/1cA5vF7M+QbwWfbgKBq3vEr7z3LffPkhsUN0REQxzvY6GHq5TPfYuoCZ0zFxgnMUJ3S5Kv+CpnxHsVS657N/4RyFYXjZzblzukvJGJbUFpl6eihbhRPes+RWtbY5bIGfGVMTlBHFwCCl9HxE/DWVVq8tVBIUi4CXU0o7zs7PUWlIeQMi4irgKoDJhwOZo1B1Y9m50t5s/rYcZeY6LHFLW9aLobDOOt010DX+hzJn1AX56UFLveUoZ26K36q2AtzZ4D6Mv7Qkvy9rzFTnzO/h7Gzb96xyjmLYVFejqqdEer8oHcXAICJGU2lVbwdeBr4NXJj799XOxusAZrVEyo5uZdOvLpSZ/gXbNmNkxiWJabIUhXtYpAcebzperT1jg4sozCneOop1pzgJ6/Ev5t/4cWMd09HKkv/SK/mOwja5qRkgMiq3Y3yt2kLNUDqKAcF5wPKU0osAEfFdKtvsqIgYWo0q+mw+eQO2ki0psUDm7mfLXgdTR9h8ppM2GPG0GKIjx48ei2taMnUBm0OekT3XtgKTizfS2wAL1eANuHDN/GzbCWNd/cv0aACKmWQjSkPySCIdC95paRXmWqAcXDRgeBY4PSKaqKSezqWSXLkHuJwK8+kjwL/ucaVtZEtT/6e8yNkyhWM6lh9pdPLIp08XaqMyopi4xKUeTHpAs3sca5Qpnfm9EW3NbnNeLEWzLuyen207psfd84YGp4G0TdxHK5ZpUmyLW1xEaaX9ba9LTVCmngYGKaUFEXEzlb17GxUp3euA7wM3RMT/rj72lT0vRnah17Y1bpERyHDxBKtwtL6XZud/IUavlF0aMnIaNTX/RGmLjZun1y/SmjTNnT7vQlLYBft25Ctuc57S4sLEx5uPzba1p/KTBZPJ1DMALuMWZW/mkdQMJetp4JBS+mPgj3d5+BnI1AzfgaFkK6sOlxv/E3J+85vFd9kOkTcCdZfMzJc7B0ASUwyF1UYUdzW6zfmda/MbNTpYpta2MPMrTD0DYI6U8DCfl3fLzbljSX4Ud+1UR0m+nJuV/VBLUKkFyoiigBgKuUzDPilUu0GnvRZxuLG5WFOIfb7dDcWZKMP9owX1ciGumc+Otnzn0HxHYeW022QdwfRpjB7hqLcndbqw7+bmy7NttWyG+Jz3THXJ/C55gBo9VGsc7DtKR1FANJJNY3yzXNoQjQClgdMmm/nMCXGZZPfYXofxSwW7p8OlBoyYHeCYabJ2c/JY1yxmnNx0OeFu9P1S+fbM/A1XU0wb8t9T27S6Ecc0S0fISZG1QOkoCohDyW50O17y4p+QqSrjWdrl0HnzZV4u6wJzpso5wuJwO6nD5b/tiVKle4SDA7I1xHbATKGzndmjO+s3XEh3iU/ON7Wd2VagcFmLK/JX2rNqgJL1VDAMB6Zl2so5DcfLE6hhPbUudaqqwzryi5+rZQPIC9OkBtKC/BtjT5RWX8tIPoxvcN18lg22ZWq+k7MNd+0jnKqAqSPZiOLZSfld5bk6bDtgh0X5hrsaOIoyoiggDgVyGXhuhorv5DbDiGS0Mq4jP69lNyGrrzOe/A13Uqfb4KY0O3ZP7hhcgPMOd44iTDc8MHJqft1JS3tL2vAs0cpvJDnApdguXDBfrb1xttv4B60zu2Q9FQu/OGQYSzLHZ069WEo122q22VhkU9x0MVXsFt6t1jabLcD5k/I33KHydc4+xaXBvsQnsm1XtTgHOlH2OphU1Rqp9fTwpHy6K8DJC/KbKBfMdieoJ8m/lgufnq/Wtoe5UsJj4FDol7yVYdkUzMZ2x11vP1OKCJooQVbKTa73Vi5Va5vhPwDKr7gU9V7IUufLSdjelYktzlGYMZ63qSHovsP95OX5jmLY7DrqK0kJD8sGtPY1QekoiodXOSR7o7O9C2POdBtF84/EF07IfYDLl5uGOPCpqofG5nfbztjsaJ2jl8iirXittqdjxvFuXOnkR/PZQD3TXDVUd7iLz9coXlJLq4FeLqvJEXLG+qCglPAoHrYxNDvfa/VyRjW4DXfONJE2cXuQsp801TGNrKMwGkgzGmXbt2y0PXZq/Zr/nmx2NOMTH81v6Bs+zTnEdbamIZhJ4+RNXy8ICo91uHtou8T1FMpaoIYRRURcCHyOiuv5ckrpml1+Pwf4e+Bk4IqUkutIrCEK7Sh6aMgOhe0MABvWzpkuHIXM3ZvCuk3f2LSGaoqTWnb2QDlL1AVu5+1qbduPcuIr+Y7CqN6Cn864eXa+FIrdnA2F2bKYrtzwTWU/csQgpZ6E6GK/y0Q0ANcC51OhYz0YEfNSSo/3MnsW+Cjwu/v+jPuGQjuKRGTnqW0uXs9paM+f/NXRIml64rtsuev3yy+zYUl1S/VQORCP2Z35juKWZlfkX2O9nGAZ283ZXssDjfnv6fkrHRtsyqT8U871r80hy0NISvqIzUJVuVaoXerpNGBpSukZgIi4gcr4hdccRUppRfV3g/BCd0bhHUXuQHb9xZcwG27HUTe5xUVdvfVe16MxaY7btBaJ1NMDTW9Ra599hJO2GCpkrUZe5u6LbkQTM0ZOftFJu28c69g9htBw/gbnKMYPzd/N17e5SEif1A0lvVZwqafWiJ2I1tdV5+lARVWo95fvOTyJf8BQaEdhYEc+WseiZjLLk7YKEiTTaLJ0FKYpznL0505wjiJ+lG9rT/E23fNQmyjyL5BzN8a6a1f1GNsTIDbnI2T946U29x0dKTvWa4b8XXNdSskNNtlPUWhHESSGqoHS+bAnSjOTeWlbfpoKoGOCSFXlszQBOFqqqpoN14r8rWmRXeLidHsSD6u1bX7d2M9YKyf5dckqvyD4dYvCN0CjSA9N5yG1tpWlH9NstZ5qIEteu9TT87CTR88b1jZIKLSjAGjIdBTD5IAWyzQxNRAzThSgY7JwFOKUDTB5tfvyHNeWXyy36rGWmTR+Qv6uNWuDo5o93uKc3E9N1kCe4hul5P2UphXZtkubJDOpMf+gYOXR7VRBqw1WM0dRm13zQeCYiGin4iCuAFnUGUAU2lEEKdsB2Pm6ajYwLgKxJ+3L2vNnIIdscrKpqqPb8iMQK2FtHcWpotfBSnK0XuBOqyoClad4m4ufKWQCrIRLY0P+XOuORY60cctM2S+i4YQe+0SNJDxSStsi4tPA7VRilK+mlB6LiD8FFqaU5kXEqcAtwGjg0oj4k5TSifv+7B6FdhRD2J7tAKzgXM92F19uHJJfcFSnT+DxlvxT34nHywE9sqfjpLflp3AmNDpRQOtAj+vILwqfOM/dlwkXuGs/TNCpl0x9k1p76otOfuadnfkHiz9u/p9qbUPt7VjoHEXXTFejKLqER0rpNuC2XR77o14/P0glJTXoKLyjGEbeCcdGFMOHuEJZj0hc2pOzyd3W21GMWJnP1Du2w/V0WOVbozt04jZ3X+zMEDNBb4E8KLSPcI6iURCZXr7A1eLquTnb3qUGK3FQC5QSHsVDJaLI29DNiQ/8h9awZCyjyjTFzZ12j1p7xDclRVs0/83qcKH+N2VEscLIqkpJKzuu9CQx6Mg2/13SdNuejXqhcWU+Fdj2C3WbvItUvbVzWgYFpaMoHkzqybb721SVcRRdW12I/eSw/JPzU42u03rG4fWT2TBCeQBf46PKXlFYT1dLa5n56WPzGT5G9RY8A6+1Id9R1HOWeOe5rs5nm0Xr3RvVL0qtp2JhCNuzT0Q2QmhlnbI36aGebe62LxsmqLeCpgsw4xjpKMRBe/QC55xHznZNcbnNluB1h078gdtAzQS9ng6301gGXsdR+bUBU/gGlx68r8GNtr1kueigBB2x1ARlRFE8BInGOrGexuBYL7m1EoDGQ91JyzghrTQ63Zlzh7B92i3dMXupsjeChrYv4sRfyJO2KGm0dbj6h9UpM1pPp652Rap/b8t3WvaeX/K0cxRjJg/SzOxycFGxUKHH5m3QNqIYIyOK0SJV1TPE3faVG/I3/2WCIQXw2DR50r5DbKCStj6X+creSFXcJ2ZaA/zKGCdQZ+qqb+d2tbR1/vMbz862veQutzmP/GD+98iQDQDEWQuAps2uN6omKCOK4sF0ZvtitkuDmJqGZY5s3ZTfWLSqxVVtF8uQ4sTDhaOQIm+nv5g/yQ9g/ti52bZWTl3LrIiaxtnSIV7Lp5R9lyAFXLJWOgrxPcqN9l+DJBx0N+ZHThXUQFuvdBTFQ5CyO7NtMdumqszQFV2EG5ovU2JVcm3vAu3fybeVqSfbVT7qckc4MHisTUZam/Md6Oi7JfX6HPc1XWOK33WcQjdFsphekn0UpkZVgTy59IXSURQT+Y7Cbfym5gBuAMwyeS2Hjcr/clpHYfV1kghAwiomyC7xCZfn5/rtgB6bX++YkO8oGueppWk6x31eNnJYvrGsURlHYVlvtr9kih7sUgNHAaSS9XTgwjoKG4ve6XAAACAASURBVDab4rfVnRrVlH9yXrvBFT7XtLjoZmHLtGzbU1tkN5+r8SpBw6NxhXKrU3Wx6HWYuN4VYSfIG2NO2stnu1niE7ryNe/PanRaT3/X8FvK3lLYwTUu9oU0BLbWYHBR0VBoR2FST7l2O2AbkcxJq55U3ec2TVFrb2xx9RIjP3LqOOko5KfxuM58R3Fcs5sBcRfnKXvDTJo4xjmKE3h8z0a9sE0Q/e+Ur/NXNucX+RsfcgeiZWe6dJ+KnGqEFLCtIbc2MujzhmqGQjsKSNlt/A2yS8Zu5kYAz47CNPIgbJMaVbKwfg/5jJpPtH9FrT1U8uKHCqrujMud5LW9L4aZNONc17syY6mz3yb6NP6Cq9Xav7JNsMFkjarnTPfZtaKTtUCKoGdo7rY5CKysOqHQjqJSV8qLFOxbZovfG0VNwxS+wXXmDhP1DPAbopklvaJZzt04SY6I/Um+6ZQL3NpHNLv3yBAUnj3FpQcn3+zqKxM68tNDtqa1TRS/h8oUjR0uZbMEtUJPw8FXpCi4o8hPPVnYmoYZoGQLq0bT6IgWt/aW7Y5psmpI/sZiNKoAOqZLR/HdfNOhsv5xcrOrrJsN1xbKJz/zPWXf2pWf2hrZ5A4Wa5rzacYTj3EptjN4QNk32saLGiARLsI/QFBoR1FP2BqFcVi269t8IWzh85FNbh7BsEPzr+WhYY5Sc0Kby8W3d4th4rKJ9yw5dOc2Ls621RtNj3MUZgrdCU3unhuWXONMt5Gf232nsl/c6Ebt1gKJUDWgAwUHjaOwkYcdsWrotL6wnp8Gs45i4RbH7jGOYrGcmX2mPFG2HyEcxQa1NHO7nQrvjY3vy7ZVCqyAVWUxg47Oxr1OEyV2yTn1598h9NGBYZcOTkQxGLWRwcZB5Cjqq13v6LSSHstL7mIEhh7inFbTsHynZaen2Xy56gGQFPoRP5KMlQvyTXVEcaYzN7Cn+L9u/L1sW8tKOv9R5ygaLx34YnGZehpARMQo4MvANCAB/xV4ErgRmAKsAN6bUtrtDmneNBsh2PDS1DTsrF/DF++Rb+moFsdFN5HTig1T1Nq2p0NtoG6kg9apOuKC/OL3Fvn+P9ie37sCcOryfFryiDucQ1x/aX6NYouMKKRqjm6KrRVq5Sgi4kLgc1SEy7+cUrpml983Av8MzKSSPH1fSmlFTZ5cYrAiis8B/55SujwihgFNwB8Cd6WUromIq4Grgd+v1RPaJrduyQYyWlLWURg2yCZ53ZZpYhhYW59rVmtbQcMftr812/b8bndatY7CFGKfkkX+b/JhZX/qof8931gwxwAaLq0j00iWHLSWVA1QqxpFRDQA1wLnA88BD0bEvJRS76LRx4GXUkodEXEF8BdAfo6zhhhwRxERhwNzoDKlJqW0FdgaEe8C5lbNvgHMp4aOwsKeGkzfhV27VVRibf3D6vEoEUH56bIqqUby4fwR0lFIDaQ5ovhtKcnzeKeyVxRW+R4Z2Qwb3W4T6Tvwh79aoJLFqMm2eRqwNKX0DEBE3AC8C3bqrnwX8NnqzzcDn4+ISCmlWlyAwWBEFO3Ai8DXIuIUYBHwG8C4lNKO6uQL0DcxPSKuAq4CGD+5fpdfT462Xds4IfvlsRGFcRRDxmxWay+XulNKwmGyWlp/M1rvyM+bjLvAhSt2HsW65sOzbce3u+LNLPLH29qa07zmi5S9pdPWApVidrZESmtE9J4MdV1K6brqzxNhpy/fc/CGk89rNimlbRHxCjAG5AyEGmAwHMVQ4M3Af0spLYiIz8HO7aEppRQRfXrN6o2+DuCEWcOzPas9xdez+G1zq2bzt9xy2yVunNy4cbJfpMu1Zo9rEuvPUkuDa+SGu/NNx1zgvueWqm0is/FvdY7ijA35jmJxi7vuL8oRsSaKqxUSqn65LqVkP3n7JQbDUTwHPJdS2iEteTMVR7EmItpSSqsjoo2MLLFpuLOOwuY/xSlDd32ba/dCaQ6maG+bFlevcI5i4wn5KZwHx8qC8ASpUyVUNkwqEeA4nE7VMjEOt6nDfRZPvDlfX2vS5TZaHfi+CI+apZ6eZ2fi85HVx/qyeS4ihgKHozuCaoMBdxQppRciYmVEHJdSehI4l0pe7nHgI8A11f//ay2f1xagbArHFKjtBmpy2raZzxSnAY4QVV7NDpE+zqw/X2hUAcySjiJEIDemR3YsN7gUi5kxYmc6nPhIvqMYc46jMVmBSt2PYmlVfaCG9NgHgWMiop2KQ7gCeP8uNvOo7IU/Bi4H7h6M+gRkOIqIOBR4B3AWlRlUW4BHge+nlB7by+f9b8C3qoynZ4CPAUOAmyLi41T0gN+bs1C9agk29WQ+PFZw0OgI2Ya79XLy23E8mW1rZ10gtYHMJncLl6m1399yvbI3chXNj7pDyPtP+Zay/zTXZtvawvpHNt+UbRuuRMVxLfmfLfCHnFodxmvhKKo1h08Dt1Ohx341pfRYRPwpsDClNA/4CvDNiFhKpWX0in1+4r3Ebh1FRPwJFScxH1hAJR10KHAscE3VifxOSulh86QppcX0nTU+16wDteM0HyiwTsjCsKSsoxgy0e0sph/B1j/WNrkC8sTZYhOSqqrjN7k6wstn5m+gOlVpbqPMV8yVXeL2kIOYX9Ifatlwl1K6jV06fFJKf9Tr518A76nJk+0j9vRW/jSl9Mf9/O5vI+IIPJ9kUGAb7mwe0hSobRrMRDf1bkKyEYtBxzj3RTYnyk3Ptaq1Vx7rqLoz3i6KFF9WS+uDcNOZ+alNLaz39nzTJDNDl+E0rexwqVogEXuR8io+drsbppS+v+tjETEEOCyl1JlSWotuTaodjHe3NQdTnAZX/LZaMab4bR2cZdSYE6i953YK3V0bxNAdSWKzyrf3j83fzc9ck88c2huYyX82Ar2347Rs25ndC/ds1AuTb3Vbye2XOudfC5QSHrtBRFwPfBLooVKEaY6Iz6WU/qqeFzeYsKcGo8f0MqPV2ubLbIXYRsvUg4lu7CY0RtLDt24SHe6j3MlZNRbiCApn9khHIfels7g323atqH8BzOPSbNuTGp1U+4i7HAOr61IpEVIDlI5i9zghpdQZER8AfkCFzroIKIyjqLcooElt2WjFjEK1BT7LknIMLLfx11OS4YgJrl9koWy8MGww2T+HTcW/vfv2bNvbG0UuCfgaH8u27eLv1dqjVzpHYdPJtUIpM94/DomIQ4B3A59PKb3aX0Pc/ooacZ8HBSaFU08dKYA1YpfzOlIu0houpvlZh/jEMy6imH3Ugj0b7UCdW7BGfDtf6K/1g+6+LNuer8e1dYg7EFnZlHr3DPWFGkp4FAq5r/hLVBRdfwbcGxFvAjrrdVG5qG/DnUtVmFSVzd3XM9S16SHDNLKNYnby2wnN+UN3tC7Q4lDmG4/Kj7Qeu8iJH564RLJ1Pp9vOkY6imFD8u/jOhkKtV8m5osAJ+CGLtUCRU49RcSxwO8Bb6LX3p9SOmdPf5vlKFJK/wD8Q68nfBZkB9Mgw7659Sx+24a7ejbz1RNWHsSm5MxGsVR0KwPwqjM3PR23ijw/wImHuxSOIabZU/l0oW1i6x8vXOYi0GO7Xd9FLVBhPclIaf/Bt4EvAv8ELm+3pz6KDwLXp5R2imWr3YHbIuJooC2lJKU5a4fciMK+ucNlTcPIkpuaA8A68tkdNkKwxW+zvmVUWRwrmv9MtzIAbc7c1G5uNxxT4Leko2gUjYujut1QrEsa8wd7WFHAlbIp8pOPfl3Z1wIFTz1tSyl9YW/+cE+veAzwUEQsolK8fpFKw10H8DYqKoZX9//nxYXvu8iPWOoZrdhTvOn6hvoOUTIjXwHaRfOfPSgMm+Yyq+uFM1/Z5Xo0Nja5SXGN0/OlKkascoOL3t2e3+vwJSnyZ535J5d/XdnXCkVNPQG3RsSngFvg9dx6SmmPg4L31EfxuYj4PHAOlXliJ1OR8HgC+FBK6dl9uepaoF5v2v7EbDCbnD3F23GVE8jPI1t5CBtp1bOYeayUk1i5NX/z37TU8f/XnOyceesFQtMof7wEABOX5Nc01l/kahQ66tt36SaNItcoqOhGQaVOsQMJOGpPf7jHI19KqQf4YfW/wsJSL+vZmW1z8fWkAdrmP7M520FEliVlMFLuKrYD/dEVp+Ybu8BJS6FMuCz/2kd/SV6MkB/pusgx8GwaVJ5xaoIiO4qUktOx6YXCJttg/+rMNt3TlsJqxqzajd+KKppiuY0oOmRntilQ241fO60XhO2Rjln+kGz+M/f9ilVONsNEIJZYMc6KPEx15rVAESU8IuKclNLdEfFLff0+pfTdPa1RaEdhYBvu7Mxsc9K2RT5Td7BCaTZ9Y5ycPXm1yc3c0Gln4eQkrDM3OPIoR3e9VY5CNWywK16RjkLskdNZrJa2n8UHp7kZIxXR631DQSOKt1EZrdUX3S4BpaPYH2GjFVPktU7oCFn8Np3floFlU2zmZGc597anw0QJNrr56bPuWkZOFvfdlT9Ul7gdVWojyvuYo+xr4SigeMXsHcKuKaX8tvpdkKv1NA74c2BCSumiiDgBeEtK6St7+8QDDVtzqOc8inrO47YnYXuKe0qI5VknZK/dvEc2lWTTZkcfle+ItBT8j91nt3uyOIjUUYD1xJ+5yOm4dmd/bfOvKXv4v9L+jUjEfkV0sYiIS4AT6TX9JaX0p3v6u9xP4NeBrwGfqf77KeBGKoM1Bg370yjUonZm29dpNlAzuwLqmzarZ4oNXJpFD9zZdUDmHmAisxcuOFyt3dqZPxtj6OfU0gyV0iaPXyRZUjVAkfsoIuKLQBOVZukvU5ma99Ocv819xa0ppZsi4g/gtelMg6PItQtyHUU9G8vAbaD1HIVq17ZOyFB1R/W4zfnxBvfFnyyihHqnC0zH8m1cUscrcQeR78kmtzOa89NJJ9/vJFxoceY26qsVBiL1FBEtVA7jU6jIJ703pfSG7siI+HfgdOA/Ukrv2MOyZ6SUTo6Ih1NKfxIRf0NF5HWPyHUUmyNiDJXCBxFxOuDGbhUMNl9u2EatUqDONMXZk7P9spkIpGGbS9+93OBO2qZA/ZI8xdv7aHStrucDam2rPmJwM5cre1NfOfkZ6Sjk6FQj7Q/wnFu+TyRC1xj3ElcDd6WUromIq6v//v0+7P6KSpSQ0934i+r/uyJiApWRWFkaBLmO4repDPo+OiLuB8aC/IQd4KhnZ/YWEQ3ZQqlRgwW3gW5tdDRCe1I7oiefTvlkgxtEZGsahplmnfOQ2W4HNdHzfzzrJNs+NPmf841thsZlwZghWVW1KGUPYI3iXcDc6s/foDKO+g2OIqV0V0TM3fXxfnBrRIyi4lz+k8rB/59y/jBXFPA/I+JtwHFAAE+mlKRsWn2Qu7nYCMGPK81f326I5lps6slIT4BjD1kGlj3FN23Od7irmt21WDqteU+tMvHbx+XPlwDZ4fxvbjff8ilRu3mbWhqOceaf4EvK/ptu+T4xgDWKcSmlHTIIL+D5aTuhOpn0rpTSy8B3IuLfgENTSlmZoVzWUwNwMZV82VDggoggpfS3e3fZtYHhNFumkR9Xmr9B23qJSfdYB2fDaFN3WNbg5LQ1S2pE/rXb12nGiYKj09p+kZNxk+JuXyNEB6UKitLM+mW3NrJv+MxF9R0p2x/EoaA1InqfOK5LKV234x8RcScwvo+/+0zvf6SU0r7O/0kpbY+Ia4EZ1X93Q/6JJdc13kolv/UI4FTE9hMYiQ3wQ3RMvtSmHkxh3cyLAJ8GM3WHdQ0uWjlWzq9Y25CfNrOCg+M6XR1pVXO+3Kzt6bDsse0rxAQgOXbafHZ/Msd1lB8n1IABRv+llB+pAWTD3bqUUr9crpRSv0PfI2JNRLSllFZHRBvYtvU+cVdE/Bfgu1UF8GzkOoojU0on++sqLupJp10rN3NTd7DUSz24qLF+ndmtXW5zfqDpLdm2tnbTICfnmhSeza3roUuGZnKKW9rMmPgmH1Jr/y5/rexH3zY4jmKAahTzqIj4XVP9/7/WYM1PUKk3b4uIX1ApI6SUUvOe/jDXUfwgIi5IKd2xDxd5QMOktqzktdF6skJ8Ng1iHJGtlwz7xZ5temNNU/6mZeirAC+3uPSgcYpW2mKR7YoT8yiOnO26oZ8UDZeLpUbVp2RD3PPPKPOaoMJ6GhCtp2uAmyLi48DPgfcCRMQs4JMppSur/76PiurVYRHxHPDxlFKfRa2U0l7ziXMdxU+AW6oFkVcRnmh/QT3VYCvr528UtrBu7G2jmKXqmo3CzqnuGjHE2YvXOqnbsZiWNTpOqonMrFTFbVys7OnI/+zOkA50oRj4vfa2yWrtnov3/47ngdJ6SimtB87t4/GFwJW9/n1W7poRcVdK6dw9PdYXcnfPvwXeAjxic1v1hOnMtsVM27FsNi2bL6+nWuWwHvc61zfkd0+fJIuw6xtdwtxIeIzY4EprK9qmKHtDpx3V6STP1zQ7CvP0Cfmbvy3a37rgPfnGS9TS2lFMPN6tzxPSvh8UTespIg6l0mvRGhGjqRz0AZqBiTlr5DqKlcCj+5OT2IF6jUK1uXszrtQ2ChlpC0sx3drg7otxWhO68occAdzb9FZlbybcWdgU3rncmW27sdmltWzd6VJuzba1kTP3C1uRAoO96LT+HWf++jl871FQradPAL8JTKAyqXQHNgKfz1kg11E8A8yPiB+w8wi9QaXHGnhRQEunzd9wzewKcAXEeqrBgqs72JrD+iYXUcxmQbbtS21uc7ZRXEfnz7Ntn2x2tGEb3RrVVq2Sa84hc93Sdizv8x932mBc6VKhfaGgWk8PADcBl6eU/jEiPgL8FyrSINfnLJD7ipdX/xtW/a9wqKcarLW3pziT1rIRxTKpD2HqDhsPdx8VW1+Z1JkfsTzS7PIU9j4OFZyAJ5uPVWvbLnGT8ruFd6u1Eczbw053o22fxN2XRzhJ2VcEsPcNAyjhUUt8CTiv6iTmAP8H+G/AdOA6MlQ2cjuz/2RfrrKe2F86sw2FsZ4nEnv6tBGF6QEw9QzwLKmhgga6vHmKWnsS9RsHb9Natu/C0IxXNrlrMXvzrCbX3b640guWjQXMVva1chQFTD01pJQ2VH9+H5XGv+9Q6dDOouDtdseKiM+nlD4dEbdSFQTsjZSSG71VY9RzFKptXHMjQt2wXyv5YGBP8WN68jehxQ2OHmn0kgClJWTSdwCzcF2/SdSbrWzKedyl7Lc15LPH7LUcfl7+zFfr4ObJSX7P3VhHtcTdoICpp4aIGJpS2kaFSXVVr99lvZg9GX0Y+DTITpj9EPUehdpKfphdTyVTrwbrnJDRV1rT7DZnS9XsPCI/BWDJDJNfdI2wy8fmd2bb+tcJGxwzaUVL/rVYOZnzhuU7Lduj89wDcuP/kTOvBQo6CvVfgB9FxDpgC3AfQER0kNmeuSdHsQwgpTQIb8mekchPEXnZ8PqNK13OFLW2kbD2o1Ddhmj0lewmNKXTCUEvbM6PWMZZBQQZ3Cwcm99fMIXlau14Wl7L7PxrsfflfdyYbbtCfs75N2dOfnBTMxTRUaSU/iwi7qIiKX5HL/bqECq1ij1iT45ibET89m4uYFBZT0G+A6inGiy4YrYNXU1ay9YcrL6O0Vey7C5TcwBYIeoOVi/JUjtN3eEs7nWLSywjn1Vl33/DNLOd2XpgRF9yegOAAtYoSCn9pI/Hsk+ge9qxGoDDeL1BY7+Ca7hzdEfLTDIRSH1rDu4Ub2oO4OoOuuYgN2fjFG039OZJrkvcpPxO6Ha5exkkKqc1U9ZiTEpu5VhZKJf6WqJJvGbYzpCBkvDYr7AnR7E6Z/D23qAqXb4QeD6l9I6IaAduAMZQaQr5UEppjwnxeo1CtadhUxQ22k12bZtia+x2RX6jr2RnOqwb64r8JoprXem6oR+a5Oi0JuobscR1iT97iuvMNj0gZ6gOOkDcxqVjZc0hq0e4Fy6s34Frdyha6qkW2JOjqGck8RtUmup36EX9BfB3KaUbqkPAPw58oVZPZtM99tRvmoXGiMI3uJOz7Sg3DBlwDtd2Zj/cNE3ZKy0pmdZaOsltciq1Jfu+bL+A6fyfujy/URCgW/ATVuVN2Xwd+WLAAMydcI+yn++W7xNFrFHUAnvaPfcoFrU3iIgjgUuAPwN+OyICOAd4f9XkG8BnqaGjsLB0WpN6GmmOZcAy1ZntipNWX8n0aTTamkOTm1yjNuc61hwAzhMSHjZzoSbWAR1Gv0nqMS1vf1O27abtjoF3+DtcdfqsCnknG/OVdd8wBJoDCbt1FL2aNGqNvwf+O7yW2B0DvFzl+UKlrNVnIBoRV1HlAbdNzo8SLD22nt2XNloxEYUV4rOnVdMlvE108QKslfO7L+7+frbtCx1uILNNVZ68Mp+Z1u36yrTT+hTX5hu7cdzcx5xs2+FDXPr2fcNuUPbKOQP/S1n3h0JKeOwzBvwVR8Q7gLUppUViKPhrqI4SvA7gxFmHptww0HYsW60fQ4+1axsG1rBu54TWNbruadOItqbZrW2ZZiNW5uf6F3Xky6ODl183jdwPT3IpNntfVDrJ+WZ1sLDSI2/HzQa369cCZepp4HAm8M6IuJhKQqAZ+Bwwqlf34JHA83tayLxpw6U8RD3HldZTiK9nqHtLbQf6uK781NbCJjdwxwoaIkQHl0pNq6MlS8oEiY/ghkXa6XwmnZSkCoaJbmwD5dmd/6Hs1zW7KLEWSIRu3jwQMOCOIqX0B8AfAFQjit9NKX0gIr5NRZzqBjJH/9Wz4c5uoOYEapviTNH25QbnhGwtxtQdrI6Q5fQj9gmb1rKnW7O86XMAn4s3/nZxi2N3mQOUUbEFGCoJWFsuct9RzWjoAwVVj91n7E+v+PeBGyLifwMPAV8Z5OtRMBGLjVZM0dYyTYz0CEASWTMbOR3d46QqXphUvxPlxOUu9dR5fP4p077/lmZsCvdWWM9Et3Nt+VjqMNoDFzgWXn8YiNRTRLQANwJTqMiBvzel9NIuNtOpEH6agR7gz1JK+a3zAoPqKFJK86mSEVJKzwCnmb/fvzqz82+l/aCN7MlPa62VQnw2xbK+Jb/Xwd7D5mdddLO4Pf+12rGsNv1trsWq5Lb+zLHkEOQxS2Ywh5aJC+Q9l20RS2Vkhmwu7AsDWKO4GrgrpXRNRFxd/ffv72LTBXw4pfR0REwAFkXE7Sklp5Gfgf0potAwndn1VIMFpwhr13bX4U6r7V2OR7+gKd+X69y63CjMJneakJ4AsHuBOZlbZpo1X/7B/KjSMqo+ytfyje9WS8MxzlwPXeImaf9GJIKe7QPiKN7F66OfvkHlQL2To+gtwZFSWhURa4GxuPFSWSi0owDTmW1rDi4lYxRh7VAcU3ew0iONkh65qik/3J8ui5lGqhvcPT+561G3+FHO3KRBLmWeW1xKW/g5DflQdQfJTeAyZ277S2qBtD3o/kV2/rU1InrnDa+rsjZzMC6ltCNX9gLsvpkqIk6jMlTO5W8zUWhHUYko8r5FNly0Rd5N4iRvBerMzOzR0gltk81fJmKx0cqyliOVvWr+kwqsS07JbyyzmLrS3RebijdifFYKfPz9oiDsMo88O9WdFJ7CUZ5rgZSCnm3Ze8m6lFK/ilQRcSd9Sxt+ZufnTCki3jAPqNc6bcA3gY+klJw+TCYK7Sj2JwaC6Y2wMhtPii+ElbCuZ69Do8w8PSl7HRSPXstmOAqrEkCUwY0VvzMy9nMso8qk+aWDs6mktU9Ndk9QCySMo9j9Uimd19/vImJNRLSllFZXHUGfvPSIaAa+D3ymL4XYWmH/2GUHADYlYyMQ0/lti7zmFG/nC9iCoE3JGdh8uZG8NrOewc9vPtkUEmRK5sGxrkHPfF507eYZYfs2t7ROmQ3GPIoUbHt1QGoU86i0CVxDP+0CETEMuAX455TSzfW8mEI7CjPo3HZm25qGUZutZ8POqA1ONmFli9ucjY5QalFL+/ndXU/kX4sslNoRoUquW+pO/VRuoCZiPXWlDG9MTeutbmk9v6LmJdscBNt7BmTbvAa4KSI+DvwceC9ARMwCPplSurL62BxgTER8tPp3H00pZc3BNii4o8invdp0j1GDBZd6sPObzSl+q9yE7Ib4zs4fZNs+0uJO5fY9ahSH+HtnK+a1ZqZNvF/ktmQqyZ601bxvWVc3jYXfH+s0RRd0yohiMHavBNQo9bTbp0lpPX2IsqaUFgJXVn/+f8D/q/vFUHBHYWDTPVYU0MyYsCkWowi7rsnVHGxKbqjIbC1tdmktPYVOyWa4foGjLXlElIaWnumK9usEmQFkSk4W+Tk93/Q+zlJLb1k62l2LG19SG2wP+MVBs22+hoPmFde7ScZIhNheB8O7t5pGepa0uI3WIb6TW921iNOtvRYzGxpQFFYbIdi06RmdP803tqoWZ+abLrShk00lOX9bO9hJfAcACu0o6tmZbemxZn17LaM68ztzVzU7CY8TcGM5zeAa2+TYvtpJLKybmn+ktDWnGS/m1z8AxfCxG6hVSR16lzCWKigPCuVbG8XZDXj4+PwBTYCcWdkPKvnugw4FdxT168y2+XLTR2HHrJra2cu48H1St9uEljblp5PsPbRaPw+11VHCw/UKgpgxsULQVwHezS3uWkwpUzYWmnkUax+X9FVZX5vQ7PjXNelEKx1F8RCk7By73UAtDXSdKArbUagbm4y+kvsUj1jr+nOWTcp3FLrmIL+ApjNXK9PKHpCfXJDvtCzr7Uypwqqu/WK3tOp1sPRVx6vQzYI1cxSv1mKhYqHQjsKgnsVpcLUBK21h1tb6ShLmWj5mdIEA2brACqF+9wGud4vLlMztvD3bdoJUMe244zl3MeJb/fBsd9Pnb5+bbyzmhQAc1uEOUO3yIOKmXfSDhNYkOxBQaEdhOrPrqQYLcmZ2t3NCaxvzq7a25tA5wTlQQwoYvcSl2JZMdbIZJnpqXSQVWN2YBh4Xf1DXbmhQ6rH3vKY7l4f1MQLUJAAAGChJREFUS/qcUNw3ZCppXJPrRNSDrmqBMvVUTORuXJYGukXOTLbFbwMjftfW7SKKpxqdbIbK9cuywGKT6EcWeSUNdPkVjhRgUpt2oI/qhgb4UL6pVmA1wU2dU0l6XG0tUDqK4sFow9sCsqWwGkXYLY2usG6ot7bmsHSSo9N2mPkV8kRpNK0AzuaefGOZBrlX9gCYwv2M5ZJRJTemdXPya1paNkPcxyETnTSxTZtaFeaaoHQUxYPpzB4lO22NYiu4D7nt+q7nF8JOxDObcxIpEPASHjO7xeQ3d8t1dKPEGO2cBnnt93B2tu3Pn5EKrGLHmDDO1WJsT49m1dUCpaMoIupXo7C8e7OZPyQ1bUyKZfMRQ9TathbTujQ/1/9Yh+vMbpTpwRELRPQke78shfX9plguBxHZa59v6g4rwi0ufLmNECyFeVAcBZSOonio3zwKi2E9+TWKlxscVddo96xodMd4Ha2I+uHiDjuWVRIYRd/Fg3OcAqs9KKg02Aa1dB+KP7uHqjvYYHV6v2MR3gDrKCwlvZ51wX6xHZ3GPBBQaEcxhJQtb2C/+Pa00tWQX/y2Tmv06vz6yu1trhv2JB5W9oY2aieQaTqtqIHcxiVqaVWLAVpvEKwqqar73baLlP3ip4Qgk9wB2o5akW1r+2hsRFHPkcL9okw9FRO5KaV6FqfB5dfrWXOwtZVLu518aOfx+XRa65w7lsp+gan5poa+CnAeRgcD1w0tlS3mi5oDgNqfXVlIpUFtY6nd+AcloigdRfFgRqHahju7mS8TA4BsSN55RP61Wx2pEU87ltS90/LTSWrqG+jc/bOX5feX2M58lUoC1w39Cbe0VWHFSCCdmp9KAvfZtfTV4dJRWBWCmqB0FMVDkLKpozbdYwurJqKwTXErGqZk22puuWMwqnSSfZ1WydRsoDaV2LFIRjfi4/Vgu6uXLF7lGFgGY6a4Q4txFHq+yGBECHuD0lEUC4nIPkFb1pOFmZk9ZoPrEn6gZUq2rS4IS6kKI9d9qZUNd1kzpcKqtZ5+5MyNuN5d9DsquW8syf9sAZgsa9sQ5yjMbBTrKOx31PQX1QxlRFE8VEah5n2JbFhraxomAgn5+TZ9F2dJeYjOY1xKzgjaTXxURjfuoK3SfZ/gS25x2RPHL+ebmj4HwIvridqNlTAfJfJaPpU0CBu/xXZqpFe+e0REC3AjMIVK1em9KaWXdrF5E5WZ2UOAQ4B/TCl9sR7XU2hHsZ3I3rhGy5qDbYozX7iXxjp5EKWv9Kj7FH9/muNempnZqsAL/PCDcsiywCVLZXFaCr89fE6+uN4dq/IFBPcGh0/L9yxWWM80xdlU0qDUHCwSDJA/uxq4K6V0TURcXf337+9isxp4S0qpOyIOAx6NiHkppZorgxbaUZiGu+HyGGC7hM0UOtvMpWiDcoOzFFYVschP1wNmfBqSfin9hByIp9RjeVSmkuTIz0nD8g8tJpUEjplkWUk2lTRoEcjA+LN3wWudk98A5rOLo0gp9b7BjVQii7qg4I6ifrAsqdau/M18UdNMtbbi9IvxoODptLM3iDBBSnhYCqsanSqjGzun4U5Td7CpJJmSM9GtbXIz6aRCRAgWrkbRGhG9dWauSyldl/m341JKOzRQXqAfIZeImAR8H+gAfq8e0QQU3lHUjx5rTysN4jthhhwBnNt9Z7bt85Pcxm9fZ4gNd9057ihsKaznkn9fLKOq82L3eXmgS3RDy/1z2JROZW/mXVgaeGGYSfWCcxTrUkr9Mi4i4k5gfB+/+sxOT5lSiog+ecwppZXAyRExAfheRNycUqq5/nqhHYXpzK53w13XiPyNxTqtESvzex3mS9kMPYVOnFfqqcAKMP5usftLdtdtDS6k2LRQOH+pqjupxRWcTf+KvecmnVSYVJJBDSU8Ukr9hqERsSYi2lJKqyOiDXafI0wprYqIR4GzgJtrc4Wvo9COwoxCXS9P8bardG1Dfs5Hd2aLXgcr1X0W97prEZucnXWg+y4MhVWmb1TNAVw39JFuactMMgN9CtHktr9hYG7BPOAjwDXV///rrgYRcSSwPqW0JSJGA28F/q4eF1N4R1GviOKkHtcmfF9D/unZbohJFFbXyiLFqUseVfam3mwL5X/Ff3fXIiis3bvyRfaAm7sud38gWmOGTXeppGNlD4hpurSppEHpXdifMHB9FNcAN0XEx4GfA+8FiIhZwCdTSldSmcP4N9W0VAB/nVKy2sRZKLSjgHzqqD0JNWxz9usa8msDc7udPMSTLfkjQnW0IvsFfnJZfmrLbion/kw2C47IN/1+kxPW23S3HM8myjE+lVS/OQ311EsqRCrJIgGvDsDTpLSePnSDU0oLgSurP/8QOLn+V1NwR7GdIdkd0fbktLXRURjNXAdTcwB4vCP/ZK47s6WEh5HN0PWP2525IUlZ9Vh76SadZFNJNg160DOT6omB66PYrzDgjqJK5/pnKnSvRIUy9rmcTsRdUenMzisM25O2pY2qoqD8bi4jf1zp5baONcGZPyKkTy/mNre47Yb+dL7pnXaogxX4nZsvrmcdhdXvMoeigz6VZFFKeAwYtgG/k1L6z4gYCSyKiB8CH2XPnYg7IRHZEYVVbDWaRiC//FLTyDT/tS9x4ye7xegCcB3rb7chgoxuHp6Z3w398wVC1wI0M8nMabCfxXoykywOyHSSQWJAJDz2Nwy4o6g2kayu/rwxIp4AJpLRifiGtUREYbXu18ii8Bk8kG27fKybU602Cnkqv3eqk80w1zL6VvmNktHNLVyWb/yYW7tPdvtuYNJsthu6ZCbtRyhTTwOPiJgCzAAWkN+JeBVwFcC4ycPqdsLZIofujF+dz+n/fpsbgqxYL1LCw8pmKMaWVWCV0c1tpn36Rbe2HGvOJDGX1dYc6slMOugjBIsy9TSwqIpYfQf4zZRSZ8TrQ9730Il4HXAdwAmzhqfcSMHSY+0pzjThmDw/wKd6/m++sewXWISTE/k//EG+8XJ3LUv/2jUY/HTBnHxjp8PIEW8WA7lxYolWoHJQJrmV6Buloxg4RMQhVJzEt1JK360+rDoRAYawPXtDXycLA0YGAWCbWN7MrgBoXpi/USyf7dJadiLeiQsEq0rOhlapJICfCVvJdrXssbLJ7SDBANFj9zcMBuspgK8AT6SU/rbXr/bYifiGtUTD3UqpBnuCTPavaM4/DRtNfwBTJ79ntpt1YGWmuUPYugybdxRGXE82WrfLcMiw6urd5Famk+qMg/D2DkZEcSbwIeCRiNck5v6QfjoRd4fKzOy8d82mnsZ1Okribc0XZNtqqQoRgCzgNLX0HDnoSM21vtYt/eNn5ro/EJ/ew05ydYE2yUwy6aQylVRg1FDrqUgYDNbTf1BpN+8LiuyeiOzObHsqGyqpmiub8+m07+NGt7hoLLO0Xt3rIG7jd8e6bmh+0N/Hoh+IdFJHk5Bqx6ceyya3gwRl6ql42M4QujLZSZaLnuRsGROxtP7Mzcxeckq+hIdNa4y+WVJYj8k3/Q5SL2mFMzcsKdslXk9mUplKKjBKemzxsJ0h2Ru07YY1+kogp9A9o5bm5lPyN1zTzwFUKkMGf5hvev0zH3Nry5B++On5tZ7jpLBeyUwq0S8OwoCw8I5iSybv8Yge1+Rk1GBB1h1kZ7ah036GP3eLy4E+PzQNel+VqSTZ5HZcc/7mb2sONgIt00kHCUp6bPGQiGx6Z9Nmd+Jb1ezahI0ibLdkA5m01skLnnKLT3bm1/OBfGNXFqio6QsYCqvVS6rnvOcylVRglMXs4mEI27Nzw92Nbu64ZUmNeDRfEfaemW9Ra6vCqqxN8zZnfsvWd+cbu1IMTHNHNZNO9BFCuZmX6ANlRFE8DGF79gawvtF1XNmNxcj33C8nv81mQb6xnEP0wp+4GaGvfFfkh1zrCkdPdnUEI65XNrmVqBkOwo9GoR1FAz3ZG/oqqTg3xepPNOebPiEnv32Mr7trEfgmH3Z/sFDYSsHWenZDN0nJz5KZVKJPlPTY4mEI2zks01E8LPWVLHuoe0a+bS6ldwcm3iLy61Lr6Xre7/7AiOtJdmyHLGoYZlK5kZeoCQaIHmvm80REM/A48L2UkpjSko9CO4ogZZ/8Xma0Wnv8UkcHuqcjv+5gqbrcLWw/4ZZefK+UbBXielZYz96Xeja5lY6lRJ8YuBrF1eTP5/lfwL31vJhCOwrTcNdo9bcdm5Lvd+RLXuspdGK//adpH3Jr/7YzNxHLTJWn8o7CpJPKSW4laoLtDNTgoqz5PBExk8pIhn9Hq6vlo9COooeG7OlvdqqYVT410hlzlvzULS4orF9DNrkZYT2A9+WbGult2JtJbtL5lyhRC+SfOVojovdp6brqmIQc7HE+T0QMAf4G+CBwXvZV7QUK7Si2MTR7trUtlK6bdpiyV9IZt6il4dJ80x8/4NRj7cjPMac+n21rIwTrKMrehRKDgvzx6OtSSv2e8iPiTvpuM/3MTk/X/3yeTwG3pZSe6z3Ppx4otKPooSG732HyateZPa8tXw0WZGf2Q2ppHvoDoQooJDaASqlM4KQh+fKxhpUEpQxGiYMLKaV+o4CIyJnP8xbgrIj4FHAYMCwiNqWUrq71tRbaUSi4xlwWynGll5kwQR5ur+XX8o2XuLW50pmbsaxeL6lMJZUoUcUe5/OklF6TSYiIjwKz6uEkoOCOwky4sykWK9c9Y4EYdCRLTl9b89F8Y8nIGHZ6p7KvZzd02btQosRr6HM+T0TMAj6ZUpJHvH1DoR3FULZla/h0tg9Ta+tNyKiwftItvf1bI/KN3dhpTmoxk4icnEiZSipx4GFgaE8ppfX0MZ8npbSQPvIAKaWvQ/06cwvtKBrpzp4zML9hrlr7bPJF/gAlnXH1pM+6tX8gbH/VLa3kQYBxou5gU0llhFBi/8fB2ZpdaEfRsC0xakOed3+4xXVm/3rPP7qLac83vbZL1BwAs98OOdON5rMDfUw6qexdKHHg4eBUBSy0o0gBPZmvYO0baci7RfO3ZdrkPfmmm/7RCRSadNKscYvU0ra/pEwnlTi4UUYUhcO2hgbWNef1OzRJ9VB+5Mx/+AUxTOF/urV5R76poukCoyQzyUhhlKmkEgceSkdROGzlkGx20km4oq2l0/41v5dvbLuhRc/l0VJYr97MpBIlDiwkBkrDY39CoR1FN4eyjI4s23O50y3uShrcMe+d+cZyTsOxJz+cbWu7octUUokSBmWNonB4lUOy50yMv0UOh5YUVj4ibKVg6ywhrmdHfpaqqiVKGJSppwMb9zvzay77TfcHi4WtVIw/lvw52GUqqUSJeqKMKAqHQ3g1n7XjWKP8Wddn9my0lxhzYb6wHrhpe2UqqUSJeqKMKAqHZjo5L7f2IEoIAJuukBTWufmmFw+5TS1tKKx27kaZSipRwqCMKAqHoRu303r3pizbf79orlv8EnkxYhaRZWDZec8lSpSoFwZuctH+hEI7CrohNyvzNT7q1s7XnAdcR7SX386PEsoIoUSJeqJMPRUP24Ff5Jne9IeGloSmx14x7sZs28mSwloWnEuU2J9Qpp6KhZHA2zJtT5Zrf9aZm3TScBm6llFCiRL7C8qIonB4afjhfHvanDzjJOZFAFx5lDI3o1ZtwblEiRL7C0pHUTj8nDfxq9u/kGmdNzJ1B949wQ22NvLbJUqUKCpK1lPh0LPkENa/dWKe8UVu7blyHkX2pD3KVFKJEsXFwLCeIqIFuJHKVPsVwHtTSi/1YdcDr+W9n00pyUaAPBTaUbD5FfjxrXm2iy5VS3eIVBKUBecSJQ4ODFjq6WrgrpTSNRFxdfXfv9+H3ZaU0vR6X0yxHQWbgf/Msrz0zZn0qCrGaM2k0lGUKHHgY8BST+/i9TbebwDz6dtRDAj2K0cRERcCnwMagC+nlK7Z/V8cBpyZtfbFfEtdS1lwLlGixBuhIorWiOit6HldSum6zL8dl1LaMaD+Beh38tqh1efYBlyTUvpe7sUZ7DeOIiIagGuB84HngAcjYl5Kqf9JPIc3w5y8YQ0n4LSbygihRIkSb4SKKNallGb198uIuBMY38evdtqsUkopIvprAX5TSun5iDgKuDsiHkkpubx5BvYbRwGcBixNKT0DEBE3UAm/+nUUozo2MHfeDVmL20luJUqUKPFG1K6YnVLq95QbEWsioi2ltDoi2oC1/azxfPX/z0TEfGAGyAJrBvYnRzERdmpZfg6YvatRRFwFXFX9Z/f34pcfzVm8LvHYwKEVWDfYFzFAOFhea/k6Bx5v2vclVt8On81VDN2X1z2PypSba6r//9ddDSJiNNCVUuqOiFYqefi/3Ifn7Bf7k6PIQjXHdx1ARCzcXWh3oOBgeZ1w8LzW8nUWEymlCwfoqa4BboqIjwM/B94LEBGzgE+mlK4Ejge+FBHbgSFUahT9p+r3AfuTo3gedhqAfWT1sRIlSpQ4qJBSWg+c28fjC4Erqz8/gFal2zsMGYgnycSDwDER0R4Rw4ArqIRfJUqUKFFiELHfRBQppW0R8Wngdir02K+mlB7bw5/lUs2KjoPldcLB81rL11miMIiU5OCFEiVKlChxUGF/Sj2VKFGiRIn9EKWjKFGiRIkSu0VhHUVEXBgRT0bE0qpo1gGJiFgREY9ExOJd5AAKjYj4akSsjYhHez3WEhE/jIinq/8fPZjXWCv081o/GxHPV9/XxRFx8WBe474iIiZFxD0R8XhEPBYRv1F9/IB8Tw82FNJR9JL7uAg4AfjliDhhcK+qrjg7pTT9QOKjA18HduWk71DMPAa4q/rvAwFf542vFeDvqu/r9JTSbQN8TbXGNuB3UkonAKcDv1b9Th6o7+lBhUI6CnrJfaSUtgI75D5KFAQppXuBDbs8/C4qSplU///uAb2oOqGf13pAIaW0OqX0n9WfNwJPUFFbOCDf04MNRXUUfcl9ZE4wKhwScEdELKrKlxzIyFXMPFDw6Yh4uJqaOmBSMhExhYrm0AIOvvf0gERRHcXBhLemlN5MJc32axGROSS82EgV3vaBzN3+AnA0MB1YDfzN4F5ObRARhwHfAX4zpdTZ+3cHwXt6wKKojuKgkfvopQ65FriFStrtQMWaqlImu1PMPBCQUlqTUupJKW0H/okD4H2NiEOoOIlvpZS+W334oHlPD2QU1VEcFHIfETEiIkbu+Bm4AMhSyy0odihmQj+KmQcKdmyeVVxGwd/XiAjgK8ATKaW/7fWrg+Y9PZBR2M7sKp3w73ld7uPPBvmSao7qMJJbqv8cClx/oLzOiPgXKqMeW4E1wB9TUYO/CZhMVTEzpVT4InA/r3UulbRTAlYAn+iVyy8cIuKtwH3AI1SGNgD8IZU6xQH3nh5sKKyjKFGiRIkSA4Oipp5KlChRosQAoXQUJUqUKFFitygdRYkSJUqU2C1KR1GiRIkSJXaL0lGUKFGiRIndonQUJQYMETGlt4Jq9bHPRsTv7uHvZkXEP1R/nhsRZ/Rj906rJBwRX4+Iy83flChxsGG/GYVaokR/qA6U3yGxPhfYBDzQh908DsDGyxIlBhtlRFFiv0FEzI+Iv4iIn0bEUxFxVvXxuRHxb1WxuU8Cv1Wd4XDWLn//0Yj4fPXnr0fEP0TEAxHxzI6oISr4fHWWyZ3AEb3+fmZE/KgqwHh7RLRFxOFV2+OqNv8SEb8yIDekRIn9BKWjKLG/YWhK6TTgN6l0ML+GlNIK4Iu8Psfhvj2s1Qa8FXgHcE31scuA46jMMfkwcAa8plP0j8DlKaWZwFeBP0spvQJ8Gvh6RFwBjE4p/dM+v8oSJQqEMvVUYiDRnwxA78d3iMktAqbs4/N9ryq693hE7JC3ngP8S0qpB1gVEXdXHz8OmAb8sCJbRAMVVVdSSj+MiPdQGZZ1yj5eU4kShUPpKEoMJNYDu85daAGW9/p3d/X/Pez757O718+xB9sAHkspveUNv4gYAhwPdFG5/uf28bpKlCgUytRTiQFDSmkTsDoizoHKPGUqI0L/QyyzERi5D5dxL/C+iGioKrieXX38SWBsRLylem2HRMSJ1d/9FpWJbe8HvlZNU5UocdCgdBQlBhofBv5nRCwG7gb+JKW0TPz9rcBlfRWzM3EL8DTwOPDPwI8BqiN1Lwf+IiJ+BiwGzqgWsa+kMg/6PiqO5n/sxfOWKFFYlOqxJUqUKFFitygjihIlSpQosVuUjqJEiRIlSuwWpaMoUaJEiRK7RekoSpQoUaLEblE6ihIlSpQosVuUjqJEiRIlSuwWpaMoUaJEiRK7xf8H/qDlscrbQsYAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 2 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 316
        },
        "id": "jHNGDVWhOhj4",
        "outputId": "01e29e40-e6ae-4886-e5de-89b6af1d1a8f"
      },
      "source": [
        "# Quality check: find and plot total energy over simulation time\n",
        "H_total = vmap(H, in_axes=(1,1,None))(x_arr, p_arr, params)\n",
        "\n",
        "plt.figure('H_over_time')\n",
        "plt.plot(sol.t, H_total)\n",
        "plt.xlabel('Time (s)')\n",
        "plt.ylabel('Energy')\n",
        "plt.ylim(0, 1.1*max(H_total))\n",
        "\n",
        "print('Maximum energy difference    = {0:.4e}'.format( (max(H_total)-min(H_total)) ))\n",
        "print('Maximum energy error         = {0:.4e}'.format( (max(H_total)-min(H_total))/jnp.mean(H_total) ))"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Maximum energy difference    = 6.1095e-07\n",
            "Maximum energy error         = 2.4438e-06\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAATR0lEQVR4nO3df7DddX3n8eerSQFbpwiScWyCJKxRG7QSvUZdttbfhlUTZwdr3DqmHXdYu2Wt63YtrjOLTaedYn+42120ZJRdti2myrptuquyDGDbmQ40N4BgwmYMqUJSLVEQl8IAgff+8f1ee7j53OQA95ubnPt8zJy55/v5fL7nvL/zhfPK98f5nFQVkiTN9kMLXYAk6fhkQEiSmgwISVKTASFJajIgJElNSxe6gPlyxhln1MqVKxe6DEk6oezcufM7VbWs1TcxAbFy5Uqmp6cXugxJOqEk+eZcfZ5ikiQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU2DBkSS9Un2JNmb5OJG/4eS7E5yW5Lrkpw10vdYklv7x/Yh65QkHW6wHwxKsgS4DHgTsB/YkWR7Ve0eGXYLMFVVDyb5BeDjwLv6voeq6tyh6pMkHdmQvyi3DthbVfsAkmwDNgI/CIiqumFk/I3AewasZ06/+me72P2331+It5akp23Nj/8Yl7z9nHl/3SFPMS0H7h5Z3t+3zeV9wJdGlk9JMp3kxiTvaK2Q5MJ+zPTBgweffsWSpB84Ln6TOsl7gCngp0eaz6qqA0nOBq5PcntV3Tm6XlVtBbYCTE1N1VN9/yGSV5JOdEMeQRwAzhxZXtG3PUGSNwIfBTZU1cMz7VV1oP+7D/gKsHbAWiVJswwZEDuA1UlWJTkJ2AQ84W6kJGuBy+nC4Z6R9tOSnNw/PwM4j5FrF5Kk4Q12iqmqDiW5CLgGWAJcUVW7kmwBpqtqO/BbwDOBzycBuKuqNgA/AVye5HG6EPvNWXc/SZIGlqqnfOr+uDI1NVXT09MLXYYknVCS7KyqqVaf36SWJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNgwZEkvVJ9iTZm+TiRv+HkuxOcluS65KcNdK3OcnX+8fmIeuUJB1usIBIsgS4DDgfWAO8O8maWcNuAaaq6ieBq4GP9+ueDlwCvBJYB1yS5LShapUkHW7II4h1wN6q2ldVjwDbgI2jA6rqhqp6sF+8EVjRP38LcG1V3VtV9wHXAusHrFWSNMuQAbEcuHtkeX/fNpf3AV96MusmuTDJdJLpgwcPPs1yJUmjjouL1EneA0wBv/Vk1quqrVU1VVVTy5YtG6Y4SVqkhgyIA8CZI8sr+rYnSPJG4KPAhqp6+MmsK0kazpABsQNYnWRVkpOATcD20QFJ1gKX04XDPSNd1wBvTnJaf3H6zX2bJOkYWTrUC1fVoSQX0X2wLwGuqKpdSbYA01W1ne6U0jOBzycBuKuqNlTVvUl+jS5kALZU1b1D1SpJOlyqaqFrmBdTU1M1PT290GVI0gklyc6qmmr1HRcXqSVJxx8DQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElS01gBkeR3kpwzdDGSpOPHuEcQdwBbk9yU5P1JTh2yKEnSwhsrIKrq01V1HvBeYCVwW5KrkrzuSOslWZ9kT5K9SS5u9L8myc1JDiW5YFbfY0lu7R/bx98kSdJ8GPsaRJIlwIv6x3eArwIfSrLtCOMvA84H1gDvTrJm1rC7gJ8Drmq8xENVdW7/2DBunZKk+bF0nEFJPgG8HbgO+I2q+uu+69Ike+ZYbR2wt6r29a+xDdgI7J4ZUFXf6Psef0rVS5IGM+4RxG3AS6vqX46Ew4x1c6yzHLh7ZHl/3zauU5JMJ7kxyTtaA5Jc2I+ZPnjw4JN4aUnS0Yx1BEF3OumFSUbb7ge+WVX3z3tVnbOq6kCSs4Hrk9xeVXeODqiqrcBWgKmpqRqoDklalMYNiE8CL6M7kgjwYmAXcGqSX6iq/9NY5wBw5sjyir5tLFV1oP+7L8lXgLXAnUdcSZI0b8Y9xfS3wNqqmqqql9N9WO8D3gR8fI51dgCrk6xKchKwCRjrbqQkpyU5uX9+BnAeI9cuJEnDGzcgXlBVu2YWqmo38KKZC9AtVXUIuAi4hu57FJ+rql1JtiTZAJDkFUn2A+8ELk8y8x4/AUwn+SpwA/Cb/XtKko6RcU8x7U7yKWDmltZ39W0nA4/OtVJVfRH44qy2/zDyfAfdqafZ6/0V8JIxa5MkDWDcI4jNwF7gg/1jH933Fx4FjvhlOUnSiemoRxD9F96+WFWvA36nMeSBea9KkrTgjnoEUVWPAY87/5IkLS7jXoN4ALg9ybXA3880VtUHBqlKkrTgxg2IL/QPSdIiMVZAVNWVSZ4BPK+q5pp7SZI0Qcb9waC3A7cCX+6Xz3UKbkmabOPe5voxukn5vgdQVbcCZw9UkyTpODBuQDzamJTPKbolaYKNe5F6V5J/DixJshr4APBXw5UlSVpo4x5B/GvgHOBh4LPA9+m+US1JmlDj3sX0IPDR/iFJWgTG/cnRFwC/DKwcXaeqXj9MWZKkhTbuNYjPA78PfBp4bLhyJEnHi3ED4lBVfWrQSiRJx5VxL1L/WZJ/leS5SU6feQxamSRpQY17BLG5//vvRtoKvywnSRNr3LuYVg1diCTp+HLEU0xJPjzy/J2z+n5jqKIkSQvvaNcgNo08/8isvvXzXIsk6ThytIDIHM9by5KkCXK0gKg5nreWJUkT5GgXqV+a5Pt0RwvP6J/TL58yaGWSpAV1xICoqiXHqhBJ0vFl3C/KSZIWGQNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqWnQgEiyPsmeJHuTXNzof02Sm5McSnLBrL7NSb7ePzbPXleSNKzBAiLJEuAy4HxgDfDuJGtmDbsL+Dngqlnrng5cArwSWAdckuS0oWqVJB1uyCOIdcDeqtpXVY8A24CNowOq6htVdRvw+Kx13wJcW1X3VtV9wLU4e6wkHVNDBsRy4O6R5f1927ytm+TCJNNJpg8ePPiUC5UkHe6EvkhdVVuraqqqppYtW7bQ5UjSRBkyIA4AZ44sr+jbhl5XkjQPhgyIHcDqJKuSnET363Tbx1z3GuDNSU7rL06/uW+TJB0jgwVEVR0CLqL7YL8D+FxV7UqyJckGgCSvSLIfeCdweZJd/br3Ar9GFzI7gC19myTpGEnVZPww3NTUVE1PTy90GZJ0Qkmys6qmWn0n9EVqSdJwDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUNGhBJ1ifZk2Rvkosb/Scn+eO+/6YkK/v2lUkeSnJr//j9IeuUJB1u6VAvnGQJcBnwJmA/sCPJ9qraPTLsfcB9VfX8JJuAS4F39X13VtW5Q9UnSTqyIY8g1gF7q2pfVT0CbAM2zhqzEbiyf3418IYkGbAmSdKYhgyI5cDdI8v7+7bmmKo6BNwPPLvvW5XkliR/nuSnWm+Q5MIk00mmDx48OL/VS9Iid7xepP4W8LyqWgt8CLgqyY/NHlRVW6tqqqqmli1bdsyLlKRJNmRAHADOHFle0bc1xyRZCpwKfLeqHq6q7wJU1U7gTuAFA9YqSZplyIDYAaxOsirJScAmYPusMduBzf3zC4Drq6qSLOsvcpPkbGA1sG/AWiVJswx2F1NVHUpyEXANsAS4oqp2JdkCTFfVduAzwB8k2QvcSxciAK8BtiR5FHgceH9V3TtUrZKkw6WqFrqGeTE1NVXT09MLXYYknVCS7KyqqVbf8XqRWpK0wAwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpKZBAyLJ+iR7kuxNcnGj/+Qkf9z335Rk5UjfR/r2PUneMmSdkqTDDRYQSZYAlwHnA2uAdydZM2vY+4D7qur5wCeAS/t11wCbgHOA9cAn+9eTJB0jQx5BrAP2VtW+qnoE2AZsnDVmI3Bl//xq4A1J0rdvq6qHq+pvgL3960mSjpGlA772cuDukeX9wCvnGlNVh5LcDzy7b79x1rrLZ79BkguBC/vFB5LseRr1ngF852msfyJabNu82LYX3ObF4uls81lzdQwZEIOrqq3A1vl4rSTTVTU1H691olhs27zYthfc5sViqG0e8hTTAeDMkeUVfVtzTJKlwKnAd8dcV5I0oCEDYgewOsmqJCfRXXTePmvMdmBz//wC4Pqqqr59U3+X0ypgNfDXA9YqSZplsFNM/TWFi4BrgCXAFVW1K8kWYLqqtgOfAf4gyV7gXroQoR/3OWA3cAj4xap6bKhae/NyquoEs9i2ebFtL7jNi8Ug25zuH+ySJD2R36SWJDUZEJKkpkUfEEebDmQSJDkzyQ1JdifZleSX+vbTk1yb5Ov939MWutb5lmRJkluS/K9+eVU/rcvefpqXkxa6xvmU5FlJrk7yf5PckeTVk76fk/yb/r/rryX5bJJTJm0/J7kiyT1JvjbS1tyv6fxev+23JXnZU33fRR0QY04HMgkOAf+2qtYArwJ+sd/Oi4Hrqmo1cF2/PGl+CbhjZPlS4BP99C730U33Mkn+E/DlqnoR8FK6bZ/Y/ZxkOfABYKqqXkx3Q8wmJm8//ze6aYdGzbVfz6e783M13ReJP/VU33RRBwTjTQdywquqb1XVzf3z/0f3obGcJ051ciXwjoWpcBhJVgBvBT7dLwd4Pd20LjBh25zkVOA1dHcHUlWPVNX3mPD9THc35jP671L9CPAtJmw/V9Vf0N3pOWqu/boR+O/VuRF4VpLnPpX3XewB0ZoO5LApPSZJP2PuWuAm4DlV9a2+69vAcxaorKH8R+DDwOP98rOB71XVoX550vb3KuAg8F/702qfTvKjTPB+rqoDwG8Dd9EFw/3ATiZ7P8+Ya7/O2+faYg+IRSXJM4H/AXywqr4/2td/QXFi7nlO8jbgnqraudC1HENLgZcBn6qqtcDfM+t00gTu59Po/sW8Cvhx4Ec5/FTMxBtqvy72gFg0U3ok+WG6cPijqvpC3/x3M4ee/d97Fqq+AZwHbEjyDbpTh6+nOz//rP5UBEze/t4P7K+qm/rlq+kCY5L38xuBv6mqg1X1KPAFun0/yft5xlz7dd4+1xZ7QIwzHcgJrz/3/hngjqr63ZGu0alONgN/eqxrG0pVfaSqVlTVSrr9en1V/SxwA920LjB52/xt4O4kL+yb3kA3G8HE7me6U0uvSvIj/X/nM9s8sft5xFz7dTvw3v5uplcB94+cinpSFv03qZP8U7pz1TPTgfz6Apc075L8E+Avgdv5h/Px/57uOsTngOcB3wR+pqpmXwg74SV5LfDLVfW2JGfTHVGcDtwCvKeqHl7I+uZTknPpLsqfBOwDfp7uH4ITu5+T/CrwLrq79W4B/gXdOfeJ2c9JPgu8lm5a778DLgH+hMZ+7YPyv9CdansQ+Pmqmn5K77vYA0KS1LbYTzFJkuZgQEiSmgwISVKTASFJajIgJElNBoTUkOTZSW7tH99OcqB//kCSTw70nh9M8t4j9L+t/0VG6ZjwNlfpKJJ8DHigqn57wPdYCtwMvGxkDqHZY9KPOa+qHhyqFmmGRxDSk5DktSO/LfGxJFcm+csk30zyz5J8PMntSb7cT29Ckpcn+fMkO5NcM8fMmq8Hbp4JhyQfSPf7Hbcl2QY/mG/nK8DbjsnGatEzIKSn5x/RfbhvAP4QuKGqXgI8BLy1D4n/DFxQVS8HrgBa39Y/j24W0hkXA2ur6ieB94+0TwM/Ne9bITUsPfoQSUfwpap6NMntdNO1fLlvvx1YCbwQeDFwbXeGiCV001LP9lye+MNGtwF/lORP6KZUmHEP3ayl0uAMCOnpeRigqh5P8mj9w0W9x+n+/wqwq6pefZTXeQg4ZWT5rXQ//vN24KNJXtKffjqlHysNzlNM0rD2AMuSvBq6adeTnNMYdwfw/H7MDwFnVtUNwK8ApwLP7Me9APhaY31p3hkQ0oD6n7K9ALg0yVeBW4F/3Bj6JbojBuhOQ/1hf9rqFuD3+p8OBXgd8L+HrVrqeJurdJxI8j+BD1fV1+fofw5wVVW94dhWpsXKgJCOE/0P/Tyn/4H6Vv8rgEer6tZjW5kWKwNCktTkNQhJUpMBIUlqMiAkSU0GhCSpyYCQJDX9f3UzijhHEmQ1AAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LO4kcX7qGBhQ"
      },
      "source": [
        "Note that energy drift (damping) is a side effect of the integration method used, though for this time span it is quite small."
      ]
    }
  ]
}